!! Parcas_40.f90 is a part of the PACIAE event generator.
!! Copyright (C) 2025 PACIAE Group.
!! PACIAE is licensed under the GNU GPL v2 or later, see LICENSE for details.
!! Open source: https://github.com/ArcsaberHep/PACIAE4
!! Author: Ben-Hao Sa, November 2002 - February 2025.

!> This is the program to deal with the parton cascade (partonic rescattering).

!!                                             By Ben-Hao at CIAE on 19/11/2002
!!                                  Last updated by An-Ke at UiO  on 04/02/2025


        subroutine parcas( time_par )
!!      Deals with the parton cascade (partonic rescattering).
!
!       It was written by Ben-Hao Sa on 19/11/2002.
!       Its input messages are in "PYJETS".
!       Its working block is "PYJETS".
!       Its output messages are in "PYJETS".
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa1/kjp21,non1,bp,iii,neve,nout,nosc
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_deadcone,i_LPM,i_diquark, &
                    ipad25,para1_1,para1_2
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        common/work7/reac(100),crose(100)
!       For the simulation control.
        COMMON/SA1_PY8/ i_mode, i_tune, KF_woDecay(1000), &
               KF_proj, KF_targ, win, energy_B, psno, b_min, b_max
        dimension pi00(4), pj00(4)
!       icol: current total number of collision pairs in the collision time list
!       N_old: the N before current collision
!       time_par: last collision time
!       icp: the icp-th collision pair
!       tcp: the collision time of the icp-th collision pair
!       i_type: event type. =1, elastic scattering; =2, inelastic; =3, shower.
!       pi00, pj00: four-momenta before collisions


        !Lei_debug
        ! Energy loss and Mean free path.
        common/mean_path/ path_sum(KSZJ), coll_sum(KSZJ)
        do i=1,N,1
            path_sum(i) = 0D0
            coll_sum(i) = 0D0
        end do
        !Lei_debug


!       Initialization.
        time_par = 0D0
        if( ABS( adj1(1) ) <= 1D-15 ) return
        if( N < 2 ) return
        call reset_eve

!       Creates the (initial) parton-parton collision time list.
        call ctlcre_par
        if( icol == 0 ) return


!-------------------------------------------------------------------------------
!-----------------------------   Parton Cascade   ------------------------------
!       Loops over sub-events within an event.
        n_loop = 0
        do while(.true.)
            N_old = N

!       Finds an event (icp) with the minimum time (tcp).
            call find_par( icp, tcp )
            time_par = tcp
            ic     = lc(1,icp)
            jc     = lc(2,icp)
            i_type = lc(6,icp)
            do i=1,4,1
                pi00(i) = P(ic,i)
                pj00(i) = P(jc,i)
            end do

!       Handles a collision (elastic or inelastic) event.
            if( i_type == 1 .OR. i_type == 2 )then
                call collis( icp )

!       Handles a final-state (medium-induced) shower event.
            else if( i_type == 3 )then
                call perform_shower( icp )
            end if

!       Performs the classical Newton motion in the Lab. frame.
            call his_p( ic, jc, ic_sh, jc_sh, pi00, pj00, tcp, i_type )

!       Updates the time list.
            if( i_type == 1 .OR. i_type == 2 )then
                call update_ctl( icp )
            else if( i_type == 3 )then
                call update_time_list_shower( icp )
            end if

            n_loop = n_loop + 1
            if( icol == 0 ) exit
        end do
!-----------------------------   Parton Cascade   ------------------------------
!-------------------------------------------------------------------------------


!       Cross sections statistics.
        do i=1,100,1
            reaci = reac(i)
            if( reaci > 0D0 ) crose(i) = crose(i) / reaci
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine reset_eve
!!      Initializes the collision time list, etc.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        LOGICAL IS_EXIST, INTO_PRS
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/shower_2/ NSH2, NSH2_max, IPTR_SH2(KSZJ), &
                         KSH2(KSZJ,8), PSH2(KSZJ,7), VSH2(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/syspar_p/rsig1,pio,tcut
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        common/work7/reac(100),crose(100)
!       For the simulation control.
        COMMON/SA1_PY8/ i_mode, i_tune, KF_woDecay(1000), &
               KF_proj, KF_targ, win, energy_B, psno, b_min, b_max
!       ic, jc: current line numbers in "PYJETS" of the collision pair
!       icol: current total number of collision pairs in the collision time list
!       reac and crose: arrays to account for the number and
!        the value of cross sections for corresponding partonic processes
!       taup(i) : formation time of particle i
!       ishp(i)=1, i-th particle will participate in the cascade
!              =0, not


        pio = 3.141592653589793D0
        ic  = 0
        jc  = 0
        ic_sh = 0
        jc_sh = 0
        N_old  = N
        NSH2 = 0
        NSH2_max = INT( KSZJ / 4 )
        ISH2 = 0
        icol  = 0
        reac  = 0D0
        crose = 0D0
        do i=1,N,1
            taup(i) = V(i,4)
            ! Excludes non-partons, historical entries and the
            !  partons that are outside of the simulation time.
            ishp(i) = 0
            KS = K(i,1)
            KF = K(i,2)
            if( .NOT.IS_EXIST(KS,i_mode) .OR. .NOT.INTO_PRS(KF) &
                .OR. V(i,4) > adj1(28) ) cycle
            ishp(i) = 1
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        LOGICAL FUNCTION INTO_PRS( KF )
!!      Determines whether a particle participates in the partonic rescattering.
        IMPLICIT NONE
        LOGICAL IS_PARTON, IS_DIQUARK
        INTEGER, INTENT(IN) :: KF


        INTO_PRS = IS_PARTON( KF ) .AND. .NOT.IS_DIQUARK( KF )


        RETURN
        END



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine ctlcre_par
!!      Creates the initial collision list.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


        ! Time resolution.
        dddt = adj1(19)
        ! The collision pair counter.
        icol = 1
        do i = 1, N-1, 1
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = i+1, N, 1
                if( ishp(j) == 0 ) cycle
                ! Forbides scatterings between q/qbar broken from a diquark.
                if( K(i,1) == 543 .AND. K(j,1) == 543 &
                    .AND. K(i,3) == K(j,3) ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Throws away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = n_jump_out + 1
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time finding.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with 1 not 0.
        icol = icol - 1


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine coij_p( ic, jc, i_fail, l_coll, t_coll, sig_t_theta_p )
!!      Calculates the collision time and kinematics of two particles.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/syspar_p/rsig1,pio,tcut
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: ic, jc
        integer :: l_coll(6)
        real(kind=8) :: t_coll
        real(kind=8) :: sig_t_theta_p(4)
        dimension ri(4), rj(4), ric(4), rjc(4)
        dimension pi(4), pj(4), pic(4), pjc(4), pij(4)
        dimension dr(3), db(3), vi(3), vj(3)
        dimension b(3)


        i_fail = 1
        l_coll = 0
        t_coll = 0D0
        sig_t_theta_p = 0D0
        t_max  = adj1(28)

        pi(4)  = P(ic,4)
        pj(4)  = P(jc,4)
        do k1=1,3,1
            pi(k1) = P(ic,k1)
            pj(k1) = P(jc,k1)
            b(k1)  = ( pi(k1) + pj(k1) ) / ( pi(4) + pj(4) )
        end do
        pic = pi
        pjc = pj
        pij = pi + pj
!       Invariant mass squared.
        eiej2 = pij(4)*pij(4) - pij(1)*pij(1) - pij(2)*pij(2) - pij(3)*pij(3)
!       Inserts the energy cut.
        dm1  = P(ic,5)
        dm2  = P(jc,5)
        ecut = SQRT(eiej2) - dm1 - dm2
        if( ecut <= 0D0 ) return
        do n1=1,4,1
            ri(n1)  = V(ic,n1)
            rj(n1)  = V(jc,n1)
        end do
        ric = ri
        rjc = rj
        KF1 = K(ic,2)
        KF2 = K(jc,2)
        call fsig( KF1, KF2, KF3, KF4, dm1, dm2, dm3, dm4, eiej2, &
                   sig, tsmp, cos_theta_s, i_proc )
        if( sig <= 0D0 ) return
        rsig1 = SQRT( sig / pio )

!       Boosts into the CM frame.
        ! Momentum.
        call lorntz( 0, b, pic, pjc )
        ! Position.
        call lorntz( 0, b, ric, rjc )

!       Finds the minimum distance and the time.
        rb = 0D0
        bb = 0D0
        rtai = 0D0
        do k1=1,3,1
            vi(k1) = pic(k1) / pic(4)
            vj(k1) = pjc(k1) / pjc(4)
        end do
        do k1=1,3,1
            dr(k1) = ric(k1) - rjc(k1) - ( vi(k1)*ric(4) - vj(k1)*rjc(4) )
            db(k1) = vi(k1)  - vj(k1)
            rb = rb + dr(k1)*db(k1)
            bb = db(k1)*db(k1) + bb
        end do
        if( bb <= 1D-10 ) return
        tcol = 0D0 - rb / bb
!       The collision should happen in the future of the CM frame.
        if( tcol < MAX( ric(4), rjc(4) ) ) return
        do ik=1,3,1
            dr(ik) = dr(ik) + tcol * db(ik)
            rtai = rtai + dr(ik) * dr(ik)
        end do
        sg = rtai
        dmin = SQRT(sg)
        if( dmin > rsig1 ) return
!       Moves along the Newton trajectory in CMS.
        do ik=1,3,1
            ric(ik) = ric(ik) + vi(ik) * ( tcol - ric(4) )
            rjc(ik) = rjc(ik) + vj(ik) * ( tcol - rjc(4) )
        end do
        ric(4) = tcol
        rjc(4) = tcol

!       Transforms back to the Lab frame (causality violation).
        call lorntz( 1, b, ric, rjc )
!       Chooses the min one for the causality violation.
        t_coll = MIN( ric(4), rjc(4) )
!       The collision should happen in the future of the Lab frame.
        if( t_coll < MAX( ri(4), rj(4), 0D0 ) ) return
!       Max time constraint (default 10000 fm/c).
        if( t_coll > t_max ) return

!       Kinematics of inelastic scatterings.
        i_proc1 = MOD( i_proc, 10 )
        if( i_proc1 == 2 .OR. i_proc1 == 6 .OR. i_proc1 ==7 )then
            i_type = 2
            pp2 = ( eiej2 - (dm3 + dm4)*(dm3 + dm4) ) &
                * ( eiej2 - (dm3 - dm4)*(dm3 - dm4) ) &
                / 4D0 / eiej2
!       Kinematics of elastic scatterings.
        else
            i_type = 1
            pp2 = pic(1)*pic(1) + pic(2)*pic(2) + pic(3)*pic(3)
        end if
!       Too small/zero momentum.
        if( pp2 < 1D-20 ) return

!       Succeeded.
        l_coll(1) = ic
        l_coll(2) = jc
        l_coll(3) = KF3
        l_coll(4) = KF4
        l_coll(5) = i_proc
        l_coll(6) = i_type
        sig_t_theta_p(1) = sig
        sig_t_theta_p(2) = tsmp
        sig_t_theta_p(3) = cos_theta_s
        sig_t_theta_p(4) = SQRT( pp2 )

        i_fail = 0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine fsig( kf1, kf2, kf3, kf4, m1, m2, m3, m4, s, &
                         sig, tsmp, cos_theta_s, i_proc  )
!!      Calculates the total cross section, decides the type of reaction,
!!       samples the t value and the scattering angle cos(theta_s).
!
!       The masses of the heavy c and b are retained (massive) in any case, and
!        other masses of light quarks might treated to zero (zero mass
!        approximation in kinematics treatment, massless).
!
!       kf1, kf2: kf codes of the incoming pair.
!       kf3, kf4: kf codes of the outgoing pair.
!       m1, m2: masses of the incoming pair.
!       m3, m4: masses of the outgoing pair.
!       s: the invariant mass squared of the colliding pair.
!       sig: the total cross section of parton kf1 colliding with kf2.
!       tsmp: the momentum transfer value sampled.
!       cos_theta_s: the scattering angle COS(theta_s).
!       i_proc: the internal processes number.
!
!       i_proc:
!        In the following internal order, we do not distinguish between q and Q.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!          xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!          xx: hQ + hQbar -> hOnia + g        (regeneration)
!         +10: with one heavy quark
!         +20: with two heavy quark
!
!       References:
!        Massless:
!            B.L. Combridge et al., Phys. Lett. B 70 (1977) 234;
!            Bin Zhang, Comput.Phys.Commun. 109 (1998) 193-206;
!            Jussi Auvinen et al., Phys.Rev.C 82 (2010) 024906;
!            Ben-Hao Sa, Comput.Phys.Commun. 183 (2012) 333–346.
!        Massive:
!            B.L. Combridge Nucl. Phys. B 151 (1979) 429-456;
!            Hamza Berrehrah et al., Phys.Rev.C 89 (2014) 5, 054901.
!            V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
!        t sampling:
!            Klaus Geiger et al., Nucl.Phys.B 369 (1992) 600-654;
!            Ben-Hao Sa et al., Comput.Phys.Commun. 183 (2012) 333–346.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_deadcone,i_LPM,i_diquark, &
                    ipad25,para1_1,para1_2
        common/sa33/smadel,ecce,secce,parecc,iparres
        integer, intent(in) :: kf1, kf2
        real(kind=8), intent(in) :: m1, m2, s
        real(kind=8) :: m3, m4
!       Note the following logical type.
        logical :: l_gg
        logical :: l_qg, l_qbarg
        logical :: l_q1q1, l_q1barq1bar
        logical :: l_q1q1bar
        logical :: l_q1q2, l_q1q2bar, l_q1barq2bar
        logical :: l_ghOnia
        real(kind=8) :: sigma_relative( 0:7 )
        real(kind=8) :: t_0_relative( 0:7 ),   t_pi_relative( 0:7 ), &
                        t_max_relative( 0:7 ), t_min_relative( 0:7 )
        real(kind=8) :: alpha_s, Lambda_QCD2, K_factor
        integer :: KF_outgoing( 2, 0:7 )


        sig  = 0D0
        tsmp = 0D0
        cos_theta_s = 0D0
        i_proc = 0
        kf3  = 0
        kf4  = 0
        dm3  = 0D0
        dm4  = 0D0
        m3   = 0D0
        m4   = 0D0
        if( s <= 0D0 ) return


!-------------------------------------------------------------------------------
!---------------------------   Process Classifying   ---------------------------
!       Classifies the incident partons.
!       kf: 1 ~ 5 -- d, u, s, c, b; 21 -- g
!           443 -- J/psi, 100443 -- Psi';
!           553 -- Upsilon, 100553 -- Upsilon';
        ikf1 = ABS(kf1)
        ikf2 = ABS(kf2)

        !Lei_debug
        ! Ignores hQ + hQ scatterings.
        if( ikf1 >= 4 .AND. ikf1 <= 5 .AND. ikf2 >= 4 .AND. ikf2 <= 5 ) return
        !Lei_debug

        l_gg = .false.
        l_qg = .false.
        l_qbarg = .false.
        l_q1q1  = .false.
        l_q1barq1bar = .false.
        l_q1q1bar    = .false.
        l_q1q2       = .false.
        l_q1q2bar    = .false.
        l_q1barq2bar = .false.
        l_ghOnia     = .false.
        ! g + g ->
        if( kf1 == 21 .AND. kf2 == 21 )then
            l_gg = .true.
        ! X + g -> X + g / X
        else if( kf1 == 21 .OR. kf2 == 21 )then
            ! q + g -> q + g
            if(      ( kf1 >= 1 .AND. kf1 <= 8 ) &
                .OR. ( kf2 >= 1 .AND. kf2 <= 8 ) )then
                    l_qg = .true.
            ! qbar + g -> qbar + g
            else if( ( kf1 <= -1 .AND. kf1 >= -8 ) &
                .OR. ( kf2 <= -1 .AND. kf2 >= -8 ) )then
                    l_qbarg = .true.
            ! g + hOnia -> hQ + hQbar
            ! Only J/psi, Psi', Upsilon and and Upsilon'.
            else if( ( ikf1 == 443 .OR. ikf1 == 100443   &
                .OR.   ikf1 == 553 .OR. ikf1 == 100553 ) &
                .OR. ( ikf2 == 443 .OR. ikf2 == 100443   &
                .OR.   ikf2 == 553 .OR. ikf2 == 100553 ) )then
                    l_ghOnia = .true.
            end if
        ! q + q ->
        else if( ikf1 >= 1 .AND. ikf1 <= 8 .AND. ikf2 >= 1 .AND. ikf2 <= 8 )then
            ! q1 + q1 -> q1 + q1
            if( kf1 == kf2 .AND. kf1 > 0 )then
                l_q1q1 = .true.
            ! q1bar + q1bar -> q1bar + q1bar
            else if( kf1 == kf2 .AND. kf1 < 0 )then
                l_q1barq1bar = .true.
            ! q1 + q1bar ->
            else if( ikf1 == ikf2 .AND. kf1*kf2 < 0 )then
                l_q1q1bar = .true.
            ! q1 + q2 -> q1 + q2
            else if( ikf1 /= ikf2 .AND. kf1 > 0 .AND. kf2 > 0 )then
                l_q1q2 = .true.
            ! q1 + q2bar -> q1 + q2bar
            else if( ikf1 /= ikf2 .AND. kf1*kf2 < 0 )then
                l_q1q2bar = .true.
            ! q1bar + q2bar -> q1bar + q2bar
            else if( ikf1 /= ikf2 .AND. kf1 < 0 .AND. kf2 < 0 )then
                l_q1barq2bar = .true.
            end if
        else
            write(*,*) "Warning! fsig: wrong KFs ! kf1, kf2 =", kf1, kf2
            return
        end if
!---------------------------   Process Classifying   ---------------------------
!-------------------------------------------------------------------------------


!-------------------------------------------------------------------------------
!--------------------------   Variable Initializing   --------------------------
!       Relative cross sections among different outgoing partons in the
!        inelastic processes.
        sigma_relative = 0D0
        t_0_relative   = 0D0
        t_pi_relative  = 0D0
        t_max_relative = 0D0
        t_min_relative = 0D0
        KF_outgoing    = 0
!       Switch of inelastic processes.
        i_inelastic = iparres
!       Selection of inelastic processes when i_inelastic = 1.
!       i_inel_proc: = 1, with only light quarks related inelastic processes
!                    = 2, with both light and heavy quarks related inel. proc.
        n_flavor_sample = 0
        if( i_inelastic == 1 )then
            if( i_inel_proc == 1 )then
                n_flavor_sample = 3
            else if( i_inel_proc == 2 )then
                n_flavor_sample = 5
            end if
        end if
!       Selection of the form of differential cross section in the LO-pQCD.
!       i_LOpQCD:
!        =1, small angle and zero-mass approximations for light quarks in
!            elastic parton-parton cross sections with the regulator.
!        =2, full LO-pQCD and massive forms for parton-parton cross sections
!            with the regulator.
!        =3, small angle and zero-mass approximations for light quarks in
!            elastic parton-parton cross sections
!            with cutoffs in the integral limits.
!        =4, full LO-pQCD and massive forms for parton-parton cross sections
!            with cutoffs in the integral limits.
!       The inelastic processes always utilize non-zero-mass full LO-QCD forms.
        i_LOpQCD = INT( adj1(20) )
        dm1 = m1
        dm2 = m2
        ! Zero-mass approximation for light quarks.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 3 )then
            if( ikf1 < 4 ) dm1 = 0D0
            if( ikf2 < 4 ) dm2 = 0D0
        end if
!       Phenomenological high-order correction factor K.
        K_factor = adj1(1)
!       Invariant mass of the colliding pair.
        sqrt_s = SQRT( s )
        ! Basic threshold energy.
        if( sqrt_s < ( m1 + m2 ) ) return
!       Effective coupling constant.
        ! Constant alpha_s.
        alpha_s  = adj1(2)
        alpha_s0 = adj1(2)
        ! Renormalization scale mu_R = Q.
        Q2 = s - dm1*dm1 - dm2*dm2
        ! Number of active quarks.
        Nf = 3
        if( Q2 > amass(4)*amass(4) ) Nf = 4
        if( Q2 > amass(5)*amass(5) ) Nf = 5
        ! QCD scale.
        Lambda_QCD2 = adj1(25)*adj1(25)
        ! Running alpha_s.
        if( alpha_s0 < 1D-10 ) alpha_s = func_alpha_s( Nf, Q2, Lambda_QCD2 )
        if( alpha_s  < 1D-10 ) return
!       t_cut: the cutoff used to regulate the 't'/'u' divergence.
        ! Constant cutoff (screening masses).
        t_cut  = adj1(3)
        t_cut0 = adj1(3)
        ! Gluon screening mass.
        dmD2   = t_cut
        ! Quark medium mass.
        dmq2   = t_cut
        ! Dynamic screening masses.
        if( t_cut0 < 1D-10 )then
            T_effictive = ABS( adj1(3) )
            ! Boltzmann or Bose-Einstein & Fermi-Dirac statistics in screening.
            ! i_statistic = 0
            i_statistic = 1
            ! Gluon screening mass.
            dmD2  = func_mD2( alpha_s, Nf, T_effictive, i_statistic )
            t_cut = dmD2
            ! Quark medium mass.
            dmq2  = func_mq2( alpha_s, T_effictive, i_statistic )
        end if
!       Momentum transfer limits (scattering angles = 0 & pi) of massless proc.
        t_0  = 0D0
        t_pi = -s
!       Note that the integral variable t is negative.
!--------------------------   Variable Initializing   --------------------------
!-------------------------------------------------------------------------------


!       Calculates the total cross section for each colliding pair.


!-------------------------------------------------------------------------------
!------------------------   Cross Section Calculating   ------------------------
!       g + g ->
        if( l_gg )then
!       1. Elastic process: g + g -> g + g.
            i_proc = 1
            dm3    = 0D0
            dm4    = 0D0
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi_relative(0) = t_min
            t_0_relative(0)  = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            t_min_relative(0) = t_min
            t_max_relative(0) = t_max
            sig_out = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                                 alpha_s, K_factor, i_proc, i_LOpQCD )
            sigma_relative(0)   = sig_out
            KF_outgoing( 1, 0 ) = 21
            KF_outgoing( 2, 0 ) = 21
!       2. Inelastic process: g + g -> q + qbar.
            if( i_inelastic == 1 )then
                i_proc = 2
                do i_flavor = 1, n_flavor_sample, 1
                    ! KF = 1 : d; KF = 2 : u.
                    i_flavor_out = i_flavor
                    if( i_flavor == 1 ) i_flavor_out = 2
                    if( i_flavor == 2 ) i_flavor_out = 1
                    ! Threshold energy.
                    if( sqrt_s >= 2D0*amass( i_flavor_out ) )then
                        dm3    = amass( i_flavor_out )
                        dm4    = amass( i_flavor_out )
                        t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                        t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                        t_pi_relative( i_flavor )  = t_min
                        t_0_relative(  i_flavor )  = t_max
                        t_max = t_max - t_cut
                        t_min = t_min + t_cut
                        t_min_relative( i_flavor ) = t_min
                        t_max_relative( i_flavor ) = t_max
                        sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                             dm3, dm4, alpha_s, K_factor, &
                                             i_proc, i_LOpQCD )
                        sigma_relative( i_flavor ) = sig_out
                        KF_outgoing( 1, i_flavor ) =  i_flavor_out
                        KF_outgoing( 2, i_flavor ) = -i_flavor_out
                    end if
                end do
            end if
!           Calculates the total corss section and determines the process to be
!            happened according to the relative cross section.
            sig = 0D0
            do i = 0, n_flavor_sample, 1
                sig = sig + sigma_relative(i)
            end do
            rand_num = PYR(1)
            prob_low = 0D0
            prob_upp = 0D0
            do i = 0, n_flavor_sample, 1
                prob_upp = prob_upp + sigma_relative(i) / sig
                if( rand_num > prob_low .AND. rand_num < prob_upp ) exit
                prob_low = prob_upp
            end do
            i_proc = 1
            if( i > 0 ) i_proc = 2
            kf3 = KF_outgoing( 1, i )
            kf4 = KF_outgoing( 2, i )
            dm3 = amass( kf3 )
            dm4 = amass( kf4 )
            t_min = t_min_relative( i )
            t_max = t_max_relative( i )
            t_pi  = t_pi_relative( i )
            t_0   = t_0_relative( i )
!       3. Elastic process: q + g -> q + g or qbar + g -> qbar + g.
        else if( l_qg .OR. l_qbarg )then
            i_proc = 3
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi   = t_min
            t_0    = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
!       4. Elastic process: q1 + q1 -> q1 + q1, q1bar + q1bar -> q1bar + q1bar.
        else if( l_q1q1 .OR. l_q1barq1bar )then
            i_proc = 4
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi   = t_min
            t_0    = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
                if( i_LOpQCD == 4 ) t_min = t_min + t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
!       q1 + q1bar ->
        else if( l_q1q1bar )then
!       5. Elastic process: q1 + q1bar -> q1 + q1bar.
            i_proc = 5
            dm3    = dm1
            dm4    = dm2
            t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi_relative(0) = t_min
            t_0_relative(0)  = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
            end if
            t_min_relative(0) = t_min
            t_max_relative(0) = t_max
            sig_out = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                                 alpha_s, K_factor, i_proc, i_LOpQCD )
            sigma_relative(0)   = sig_out
            KF_outgoing( 1, 0 ) = kf1
            KF_outgoing( 2, 0 ) = kf2
            if( i_inelastic == 1 )then
                ! Non-zero mass.
                dm1 = m1
                dm2 = m2
                ! Running alpha_s.
                if( alpha_s0 < 1D-10 )then
                    Q2 = s - dm1*dm1 - dm2*dm2
                    Nf = 3
                    if( Q2 > amass(4)*amass(4) ) Nf = 4
                    if( Q2 > amass(5)*amass(5) ) Nf = 5
                    alpha_s = func_alpha_s( Nf, Q2, Lambda_QCD2 )
                end if
!       6. Inelastic process: q1 + q1bar -> q2 + q2bar.
                i_proc = 6
                do i_flavor = 1, n_flavor_sample, 1
                    ! KF = 1 : d; KF = 2 : u.
                    i_flavor_out = i_flavor
                    if( i_flavor == 1 ) i_flavor_out = 2
                    if( i_flavor == 2 ) i_flavor_out = 1
                    ! Threshold energy.
                    if( sqrt_s >= 2D0*amass( i_flavor_out ) )then
                        if( ikf1 == i_flavor_out ) cycle
                        dm3    = amass( i_flavor_out )
                        dm4    = amass( i_flavor_out )
                        t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                        t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                        t_pi_relative( i_flavor )  = t_min
                        t_0_relative(  i_flavor )  = t_max
                        t_min_relative( i_flavor ) = t_min
                        t_max_relative( i_flavor ) = t_max
                        sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                             dm1, dm3, alpha_s, K_factor, &
                                             i_proc, i_LOpQCD )
                        sigma_relative( i_flavor ) = sig_out
                        KF_outgoing( 1, i_flavor ) = i_flavor_out * SIGN(1,kf1)
                        KF_outgoing( 2, i_flavor ) = i_flavor_out * SIGN(1,kf2)
                    end if
                end do
!       7. Inelastic process: q1 + q1bar -> g + g.
                i_proc = 7
                i_flavor = 7
                n_flavor_sample = 7
                dm3    = 0D0
                dm4    = 0D0
                t_min  = func_t_pi( s, dm1, dm2, dm3, dm4 )
                t_max  = func_t_0(  s, dm1, dm2, dm3, dm4 )
                t_pi_relative( i_flavor )  = t_min
                t_0_relative(  i_flavor )  = t_max
                t_max = t_max - t_cut
                t_min = t_min + t_cut
                t_min_relative( i_flavor ) = t_min
                t_max_relative( i_flavor ) = t_max
                sig_out = get_sigma( t_min, t_max, t_cut, s,      &
                                     dm1, dm2, alpha_s, K_factor, &
                                     i_proc, i_LOpQCD )
                sigma_relative( i_flavor ) = sig_out
                KF_outgoing( 1, i_flavor ) = 21
                KF_outgoing( 2, i_flavor ) = 21
            end if
!           Calculates the total corss section and determines the process to be
!            happened according to the relative cross section.
            sig = 0D0
            do i = 0, n_flavor_sample, 1
                sig = sig + sigma_relative(i)
            end do
            rand_num = PYR(1)
            prob_low = 0D0
            prob_upp = 0D0
            do i = 0, n_flavor_sample, 1
                prob_upp = prob_upp + sigma_relative(i) / sig
                if( rand_num > prob_low .AND. rand_num < prob_upp ) exit
                prob_low = prob_upp
            end do
            i_proc = 5
            if( i == 7 )then
                i_proc = 7
            else if( i > 0 )then
                i_proc = 6
            end if
            kf3 = KF_outgoing( 1, i )
            kf4 = KF_outgoing( 2, i )
            dm3 = amass( kf3 )
            dm4 = amass( kf4 )
            t_min = t_min_relative( i )
            t_max = t_max_relative( i )
            t_pi  = t_pi_relative( i )
            t_0   = t_0_relative( i )
!       8. Elastic process: q1 + q2 -> q1 + q2 / q1    + q2bar -> q1    + q2bar
!                                              / q1bar + q2bar -> q1bar + q2bar
        else if( l_q1q2 .OR. l_q1q2bar .OR. l_q1barq2bar )then
            i_proc = 8
            dm3   = dm1
            dm4   = dm2
            t_min = func_t_pi( s, dm1, dm2, dm3, dm4 )
            t_max = func_t_0(  s, dm1, dm2, dm3, dm4 )
            t_pi  = t_min
            t_0   = t_max
            if( i_LOpQCD > 2 )then
                t_max = t_max - t_cut
            end if
            sig = get_sigma( t_min, t_max, t_cut, s, dm3, dm4, &
                             alpha_s, K_factor, i_proc, i_LOpQCD )
            kf3 = kf1
            kf4 = kf2
        end if
!------------------------   Cross Section Calculating   ------------------------
!-------------------------------------------------------------------------------


        if( sig <= 0D0 ) return
!       0.0389379 (~0.04) is the transformation factor from GeV^-2 to fm^2.
        sig = sig * 0.0389379D0
!       tsmp: the t value sampled.
        tsmp = sample_t( t_min, t_max, t_cut, s, dm1, dm2, dm3, dm4, &
                         i_proc, i_LOpQCD )
!       Symmetric scattering angle.
        ! if( PYR(1) > 0.5D0 ) tsmp = t_min + t_max - tsmp
!       The cosine of the scattering angle theta_s.
        ! cos_theta_s = 1D0 - 2D0 * ( t_0 - tsmp ) / ( t_0 - t_pi )
!       Equivalently.
        cos_theta_s = 1D0 + 2D0*s*tsmp / ( s - dm1 - dm2 ) / ( s - dm1 + dm2 )
        m3 = dm3
        m4 = dm4

        !Lei_debug
        ! call PYFILL( i_proc, tsmp, 1D0 )
        ! Special treatment for the massive processes.
        if( dm1+dm2+dm3+dm4 > 1D-5 )then
            tsmp        = t_min
            cos_theta_s = t_max
        end if
        !Lei_debug

!       Labels heavy quarks realted proceses.
        ikf3 = ABS(kf3)
        ikf4 = ABS(kf4)
        ! Double-heavy, 20+
        if(      ( ikf1 > 3 .AND. ikf2 > 3 .AND. ikf1 /= 21 .AND. ikf2 /= 21 ) &
            .OR. ( ikf3 > 3 .AND. ikf4 > 3 .AND. ikf3 /= 21 .AND. ikf4 /= 21 ) &
            )then
            i_proc = i_proc + 20
        ! Single-heavy, 10+
        else if( ( ikf1 > 3 .AND. ikf1 /= 21 ) &
            .OR. ( ikf2 > 3 .AND. ikf2 /= 21 ) )then
            i_proc = i_proc + 10
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function amass(kf)
!!      Mass of a parton.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        common/sa18/i_deex,n_deex_step,i_pT_coal,i_pT_endpoint,a_FF,aPS_c,aPS_b
        common/sa24/adj1(40),nnstop,non24,zstop


!       Hadronization model.
        i_had_model = INT( adj1(12) )
        i_mass = 0
!       For the debug, uses the current algebra mass.
        if( i_had_model == 1 .AND. i_deex > 99 ) i_mass = 3
        amass = 0D0
        r     = 0D0
        k1    = ABS(kf)
        select case( k1 )
        case( 21 )
        case( 2 )
            r = 0.33D0
            if( i_mass == 3 ) r = 0.0056D0
        case( 1 )
            r = 0.33D0
            if( i_mass == 3 ) r = 0.0099D0
        case( 3 )
            r = 0.5D0
            if( i_mass == 3 ) r = 0.199D0
        case( 4 )
            r = 1.5D0
            if( i_mass == 3 ) r = 1.23D0
        case( 5 )
            r = 4.8D0
            if( i_mass == 3 ) r = 4.17D0
        case( 6 )
            r = 175D0
            if( i_mass == 3 ) r = 165D0
        end select
        amass = r


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function get_sigma( t_min, t_max, t_cut, s, m1, m2, &
                            alpha_s, K_factor, i_proc, i_LOpQCD )
!!      Total cross section of q / qbar /g + q / qbar / g .
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        real(kind=8), parameter :: pi=3.141592653589793D0
        real(kind=8), intent(in) :: t_min, t_max, t_cut, s, m1, m2
        real(kind=8), intent(in) :: alpha_s, K_factor
        integer, intent(in) :: i_proc, i_LOpQCD
        real(kind=8) :: sigma
        real(kind=8) :: dME_PF0
        real(kind=8) :: dME_PF1_stu, dME_PF1_sut
        real(kind=8) :: dME_PF1_ust, dME_PF1_uts
        real(kind=8) :: dME_PF2_stu, dME_PF2_ust
        real(kind=8) :: dME_PF3_stu, dME_PF3_tsu
        real(kind=8) :: dME_PF0_REG
        real(kind=8) :: dME_PF1_stu_REG, dME_PF1_sut_REG
        real(kind=8) :: dME_PF1_ust_REG, dME_PF1_uts_REG
        real(kind=8) :: dME_PF2_stu_REG, dME_PF2_ust_REG
        real(kind=8) :: dME_PF3_stu_REG, dME_PF3_tsu_REG
        real(kind=8) :: get_sigma


        get_sigma = 0D0
        sigma = 0D0
        if( t_min > t_max ) return

!       Elastic processes.

!       Small angle and zero-mass approximations for light quarks in
!        elastic parton-parton cross sections with the regulator.
        if( i_LOpQCD == 1 )then
            select case( i_proc )
            case(1)
                sigma = 9D0/2D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
            case(3)
                sigma = 2D0     * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF3_tsu_REG( s, t_max, MAX( m1, m2 ), t_cut ) &
                          - dME_PF3_tsu_REG( s, t_min, MAX( m1, m2 ), t_cut )
                    sigma = - 8D0 / 3D0 * sigma
                end if
            case( 4 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_stu_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF1_sut_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_sut_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF2_stu_REG( s, t_max, m1, t_cut ) &
                          - dME_PF2_stu_REG( s, t_min, m1, t_cut )
                    ! Symmetry factor for identical particles in final state.
                    sigma = sigma / 2D0
                end if
            case( 5 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_ust_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_ust_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF1_uts_REG( s, t_max, m1, m1, t_cut ) &
                          - dME_PF1_uts_REG( s, t_min, m1, m1, t_cut ) &
                          + dME_PF2_ust_REG( s, t_max, m1, t_cut ) &
                          - dME_PF2_ust_REG( s, t_min, m1, t_cut )
                end if
            case( 8 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/(s + t_cut) )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu_REG( s, t_max, m1, m2, t_cut ) &
                          - dME_PF1_stu_REG( s, t_min, m1, m2, t_cut )
                end if
            end select

!       Full LO-pQCD and massive forms for parton-parton cross sections
!        with the regulator.
        else if( i_LOpQCD == 2 )then
            select case( i_proc )
            case(1)
                sigma = dME_PF0_REG( s, t_max, t_cut ) &
                      - dME_PF0_REG( s, t_min, t_cut )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(3)
                sigma = dME_PF3_tsu_REG( s, t_max, MAX( m1, m2 ), t_cut ) &
                      - dME_PF3_tsu_REG( s, t_min, MAX( m1, m2 ), t_cut )
                sigma = - 8D0 / 3D0 * sigma
            case(4)
                sigma = dME_PF1_stu_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_stu_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF1_sut_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_sut_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF2_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF2_stu_REG( s, t_min, m1, t_cut )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(5)
                sigma = dME_PF1_ust_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_ust_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF1_uts_REG( s, t_max, m1, m1, t_cut ) &
                      - dME_PF1_uts_REG( s, t_min, m1, m1, t_cut ) &
                      + dME_PF2_ust_REG( s, t_max, m1, t_cut ) &
                      - dME_PF2_ust_REG( s, t_min, m1, t_cut )
            case(8)
                sigma = dME_PF1_stu_REG( s, t_max, m1, m2, t_cut ) &
                      - dME_PF1_stu_REG( s, t_min, m1, m2, t_cut )
            end select

!       Small angle and zero-mass approximations for light quarks in elastic
!        parton-parton cross sections with cutoffs in the integral limits.
        else if( i_LOpQCD == 3 )then
            select case( i_proc )
            case(1)
                sigma = 9D0/2D0 * s**2 * ( 1D0/t_cut - 1D0/s )
            case(3)
                sigma = 2D0     * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF3_tsu( s, t_max, MAX( m1, m2 ) ) &
                          - dME_PF3_tsu( s, t_min, MAX( m1, m2 ) )
                    sigma = - 8D0 / 3D0 * sigma
                end if
            case( 4 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu( s, t_max, m1, m1 ) &
                          - dME_PF1_stu( s, t_min, m1, m1 ) &
                          + dME_PF1_sut( s, t_max, m1, m1 ) &
                          - dME_PF1_sut( s, t_min, m1, m1 ) &
                          + dME_PF2_stu( s, t_max, m1 ) &
                          - dME_PF2_stu( s, t_min, m1 )
                    ! Symmetry factor for identical particles in final state.
                    sigma = sigma / 2D0
                end if
            case( 5 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_ust( s, t_max, m1, m1 ) &
                          - dME_PF1_ust( s, t_min, m1, m1 ) &
                          + dME_PF1_uts( s, t_max, m1, m1 ) &
                          - dME_PF1_uts( s, t_min, m1, m1 ) &
                          + dME_PF2_ust( s, t_max, m1 ) &
                          - dME_PF2_ust( s, t_min, m1 )
                end if
            case( 8 )
                sigma = 8D0/9D0 * s**2 * ( 1D0/t_cut - 1D0/s )
                ! Non-zero mass for heavy-flavor quarks.
                if( m1+m2 > 1D-1 )then
                    sigma = dME_PF1_stu( s, t_max, m1, m2 ) &
                          - dME_PF1_stu( s, t_min, m1, m2 )
                end if
            end select

!       Full LO-pQCD and massive forms for parton-parton cross sections
!        with cutoffs in the integral limits.
        else
            select case( i_proc )
            case(1)
                sigma = dME_PF0( s, t_max ) &
                      - dME_PF0( s, t_min )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(3)
                sigma = dME_PF3_tsu( s, t_max, MAX( m1, m2 ) ) &
                      - dME_PF3_tsu( s, t_min, MAX( m1, m2 ) )
                sigma = - 8D0 / 3D0 * sigma
            case(4)
                sigma = dME_PF1_stu( s, t_max, m1, m1 ) &
                      - dME_PF1_stu( s, t_min, m1, m1 ) &
                      + dME_PF1_sut( s, t_max, m1, m1 ) &
                      - dME_PF1_sut( s, t_min, m1, m1 ) &
                      + dME_PF2_stu( s, t_max, m1 ) &
                      - dME_PF2_stu( s, t_min, m1 )
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            case(5)
                sigma = dME_PF1_ust( s, t_max, m1, m1 ) &
                      - dME_PF1_ust( s, t_min, m1, m1 ) &
                      + dME_PF1_uts( s, t_max, m1, m1 ) &
                      - dME_PF1_uts( s, t_min, m1, m1 ) &
                      + dME_PF2_ust( s, t_max, m1 ) &
                      - dME_PF2_ust( s, t_min, m1 )
            case(8)
                sigma = dME_PF1_stu( s, t_max, m1, m2 ) &
                      - dME_PF1_stu( s, t_min, m1, m2 )
            end select
        end if

!       Inelastic processes.

!       Full LO-pQCD and massive forms for inelastic processes
!        with the regulator.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 2 )then
            select case( i_proc )
            case(2)
                sigma = dME_PF3_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF3_stu_REG( s, t_min, m1, t_cut )
            case(6)
                sigma = dME_PF1_ust_REG( s, t_max, m1, m2, t_cut ) &
                      - dME_PF1_ust_REG( s, t_min, m1, m2, t_cut )
            case(7)
                sigma = dME_PF3_stu_REG( s, t_max, m1, t_cut ) &
                      - dME_PF3_stu_REG( s, t_min, m1, t_cut )
                sigma = 64D0 / 9D0 * sigma
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            end select

!       Full LO-pQCD and massive forms for inelastic processes
!        with cutoffs in the integral limits.
        else
            select case( i_proc )
            case(2)
                sigma = dME_PF3_stu( s, t_max, m1 ) &
                      - dME_PF3_stu( s, t_min, m1 )
            case(6)
                sigma = dME_PF1_ust( s, t_max, m1, m2 ) &
                      - dME_PF1_ust( s, t_min, m1, m2 )
            case(7)
                sigma = dME_PF3_stu( s, t_max, m1 ) &
                      - dME_PF3_stu( s, t_min, m1 )
                sigma = 64D0 / 9D0 * sigma
                ! Symmetry factor for identical particles in the final state.
                sigma = sigma / 2D0
            end select
        end if

        if( sigma > 0D0 )  get_sigma = K_factor * sigma * pi * alpha_s**2 &
                                     / (s - (m1 + m2)**2 ) / (s - (m1 - m2)**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function sample_t( t_min, t_max, t_cut, s, m1, m2, m3, m4, &
                           i_proc, i_LOpQCD ) result(t)
!!      Samples t from the differential cross section distribution dsigma/dt.
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        real(kind=8), intent(in) :: t_min, t_max, t_cut, s, m1, m2, m3, m4
        integer, intent(in) :: i_proc, i_LOpQCD
        real(kind=8) :: t, rand_num, rand_num_1, rand_num_2
        real(kind=8) :: dsigma_dt, dsigma_dt_max
        real(kind=8) :: dsigma_dt_func, PYR
        real(kind=8) :: prob_accept
        ! real(kind=8) :: get_sigma
        ! real(kind=8) :: weight1, weight2, max1, max2, max_seperate, t_separate
        ! real(kind=8) :: norm1, norm2, M_norm


        t = 0D0

!       Elastic processes wiht the small angle and zero-mass light quarks.
        ! Direct sampling method.
        select case( i_proc )
        case( 1, 3, 4, 5, 8 )
            select case( i_LOpQCD )
            ! Cross sections regularized by the regulator.
            case(1)
                ! With the common form of f(t) = 1/(t - t_cut)^2.
                rand_num = PYR(1)
                t = t_cut * ( 1D0 - ( s + t_cut ) / ( s*rand_num + t_cut ) )
                if( m1+m2 < 1D-1 ) return
            ! Cross sections regularized by cutoffs in the integral limits.
            case(3)
                ! With the common form of f(t) = 1/t^2.
                rand_num = PYR(1)
                t = t_min * t_max / ( t_max + (t_min - t_max)*rand_num )
                if( m1+m2 < 1D-1 ) return
            end select
        end select

        ! Just the simple sampling.
        return


!       Full LO-pQCD forms for el. and inel. cross sections. Rejection sampling.
        do while(.true.)
            rand_num_1 = PYR(1)
            t = t_min + rand_num_1 * ( t_max - t_min )
            dsigma_dt     = dsigma_dt_func( s, t, t_cut, m1, m2,    &
                                            m3, m4, i_proc, i_LOpQCD )
            dsigma_dt_max = MAX(                                    &
                        dsigma_dt_func( s, t_min, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ), &
                        dsigma_dt_func( s, t_max, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ) )
            prob_accept = dsigma_dt / dsigma_dt_max
            rand_num_2  = PYR(1)
            if( rand_num_2 <= prob_accept ) exit
        end do


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dsigma_dt_func( s, t, t_cut, m1, m2, m3, m4, &
                 i_proc, i_LOpQCD )
!!      The differential cross section f(t) = dsigma/dt.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
!       In this function, we do not distinguish between q and Q.
!       i_proc: internal order number of the process.
!           1: g     + g     -> g     + g
!           2: g     + g     -> q1    + q1bar   (inelastic)
!           3: q     + g     -> q     + g
!              qbar  + g     -> qbar  + g
!           4: q1    + q1    -> q1    + q1
!              q1bar + q1bar -> q1bar + q1bar
!           5: q1    + q1bar -> q1    + q1bar
!           6: q1    + q1bar -> q2    + q2bar   (inelastic)
!           7: q1    + q1bar -> g     + g       (inelastic)
!           8: q1    + q2    -> q1    + q2
!              q1    + q2bar -> q1    + q2bar
!              q1bar + q2bar -> q1bar + q2bar
!#TODO(Lei20241113): Heavy onia gluon dissociation.
!            xx: hOnia + g  -> hQ + hQbar       (dissociation)
!#TODO(Lei20241113): Heavy onia gluon regeneration.
!            xx: hQ + hQbar -> hOnia + g        (regeneration)
        implicit none
        integer, intent(in) :: i_proc, i_LOpQCD
        real(kind=8), intent(in) :: s, t, t_cut, m1, m2, m3, m4
        real(kind=8) :: u, dsigma_dt_func
        real(kind=8) :: dME_G0,     dME_G1,     dME_G2,     dME_G3
        real(kind=8) :: dME_G0_REG, dME_G1_REG, dME_G2_REG, dME_G3_REG


        u = m1**2 + m2**2 + m3**2 + m4**4 - s - t
        dsigma_dt_func = 0D0

!       With regulator.
        if( i_LOpQCD == 1 .OR. i_LOpQCD == 2 )then
            select case( i_proc )
            case(1)
                dsigma_dt_func = dME_G0_REG( s, t, u, t_cut )
            case(2)
                dsigma_dt_func = dME_G3_REG( s, t, u, m3, t_cut )
            case(3)
                dsigma_dt_func = - 8D0/3D0 &
                               * dME_G3_REG( t, s, u, MAX(m1,m2), -t_cut )
            case(4)
                dsigma_dt_func = dME_G1_REG( s, t, u, m1,    m1, t_cut ) &
                               + dME_G1_REG( s, u, t, m1,    m1, t_cut ) &
                               + dME_G2_REG( s, t, u, m1, t_cut, t_cut )
            case(5)
                dsigma_dt_func = dME_G1_REG( u, s, t, m1,     m1, -t_cut ) &
                               + dME_G1_REG( u, t, s, m1,     m1,  t_cut ) &
                               + dME_G2_REG( u, s, t, m1, -t_cut,  t_cut )
            case(6)
                dsigma_dt_func = dME_G1_REG( u, s, t, m1, m3, -t_cut )
            case(7)
                dsigma_dt_func = 64D0 / 9D0 * dME_G3_REG( s, t, u, m1, t_cut )
            case(8)
                dsigma_dt_func = dME_G1_REG( s, t, u, m1, m2, t_cut )
            end select
!       Without regulator.
        else
            select case( i_proc )
            case(1)
                dsigma_dt_func = dME_G0( s, t, u )
            case(2)
                dsigma_dt_func = dME_G3( s, t, u, m3 )
            case(3)
                dsigma_dt_func = - 8D0/3D0 * dME_G3( t, s, u, MAX(m1,m2) )
            case(4)
                dsigma_dt_func = dME_G1( s, t, u, m1, m1 ) &
                               + dME_G1( s, u, t, m1, m1 ) &
                               + dME_G2( s, t, u, m1 )
            case(5)
                dsigma_dt_func = dME_G1( u, s, t, m1, m1 ) &
                               + dME_G1( u, t, s, m1, m1 ) &
                               + dME_G2( u, s, t, m1 )
            case(6)
                dsigma_dt_func = dME_G1( u, s, t, m1, m3 )
            case(7)
                dsigma_dt_func = 64D0 / 9D0 * dME_G3( s, t, u, m1 )
            case(8)
                dsigma_dt_func = dME_G1( s, t, u, m1, m2 )
            end select
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_alpha_s( Nf, Q2, Lambda_QCD2 )
!!      First-order running strong couping constant alpha_s.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, intent(in) :: Nf
        real(kind=8), intent(in) :: Q2, Lambda_QCD2
        real(kind=8) :: func_alpha_s


        func_alpha_s = 0D0
        if( Q2 <= Lambda_QCD2 ) return
        func_alpha_s = 12D0 * pi / ( 33D0 - 2D0*Nf ) / LOG( Q2 / Lambda_QCD2 )
        ! Upper limit ?
        ! if( func_alpha_s > 1D0 ) func_alpha_s = 1D0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_mD2( alpha_s, Nf, T, i_statistic )
!!      Debye mass of the gluon within a thermal medium.
!       T: temperature
!       i_statistic: Bose-Einstein and Fermi-Dirac statistics or Boltzmann
!       Cf. Zhe Xu, Phys.Rev.C 71 (2005) 064901.
!           Stephen M.H. Wong, Phys.Rev.C 54 (1996) 2588-2599.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, parameter :: Nc = 3
        real(kind=8), intent(in) :: alpha_s, T
        integer, intent(in) :: Nf, i_statistic
        real(kind=8) :: func_mD2


!       Quantum statistics case, i.e. Bose-Einstein and Fermi-Dirac
!        distributions for gluons and quarks, respectively.
        if( i_statistic == 1 )then
            func_mD2 = 4D0 * pi * alpha_s / 3D0 * ( 1D0*Nc + Nf/2D0 ) * T*T

!       Classical statistics case, i.e. Boltzmann distribution for
!        both gluons and quarks.
        else
            func_mD2 = 8D0 * alpha_s / pi * ( 1D0*Nc + 1D0*Nf ) * T*T
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_mq2( alpha_s, T, i_statistic )
!!      Quark medium mass within a thermal medium.
!       T: temperature
!       i_statistic: Bose-Einstein and Fermi-Dirac statistics or Boltzmann
!       Cf. Zhe Xu, Phys.Rev.C 71 (2005) 064901.
!           Stephen M.H. Wong, Phys.Rev.C 54 (1996) 2588-2599.
        implicit none
        real(kind=8), parameter :: pi = 3.141592653589793D0
        integer, parameter :: Nc = 3
        real(kind=8), intent(in) :: alpha_s, T
        integer, intent(in) :: i_statistic
        real(kind=8) :: func_mq2


!       Quantum statistics case, i.e. Bose-Einstein and Fermi-Dirac
!        distributions for gluons and quarks, respectively.
        if( i_statistic == 1 )then
            func_mq2 = pi * alpha_s / 4D0 * ( 1D0*Nc*Nc - 1D0 ) / Nc * T*T

!       Classical statistics case, i.e. Boltzmann distribution for
!        both gluons and quarks.
        else
            func_mq2 = 2D0 * alpha_s / pi * ( 1D0*Nc*Nc - 1D0 ) / Nc * T*T
        end if


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_t_0( s, m1, m2, m3, m4 )
!!      t_pi < t < t_0: the momentum transfer interval without cutoff.
        implicit none
        real(kind=8), intent(in) :: s, m1, m2, m3, m4
        real(kind=8) :: func_t_0, func_lambda, pow2


        func_t_0 = 1D0 / 4D0 / s * pow2( ( m1*m1 - m2*m2 - m3*m3 + m4*m4 ) )  &
                 - 1D0 / 4D0 / s * pow2( ( func_lambda( s, m1*m1, m2*m2 )     &
                 -                         func_lambda( s, m3*m3, m4*m4 ) ) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_t_pi( s, m1, m2, m3, m4 )
!!      t_pi < t < t_0: the momentum transfer interval without cutoff.
        implicit none
        real(kind=8), intent(in) :: s, m1, m2, m3, m4
        real(kind=8) :: func_t_pi, func_lambda, pow2


        func_t_pi = 1D0 / 4D0 / s * pow2( ( m1*m1 - m2*m2 - m3*m3 + m4*m4 ) )  &
                  - 1D0 / 4D0 / s * pow2( ( func_lambda( s, m1*m1, m2*m2 )     &
                  +                         func_lambda( s, m3*m3, m4*m4 ) ) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function func_lambda( a, b, c )
!!      Kallen (triangle, lambda) function.
        implicit none
        real(kind=8), intent(in) :: a, b, c
        real(kind=8) :: func_lambda


        func_lambda = SQRT( a*a + b*b + c*c - 2D0*a*b - 2D0*a*c - 2D0*b*c )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G0( s, t, u )
!!      The Matrix Element component G0.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u
        real(kind=8) :: dME_G0


        dME_G0 = 9D0/2D0 * ( 3D0 - u*t/s**2 - s*u/t**2 - t*s/u**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G0_REG( s, t, u, t_cut )
!!      The Matrix Element component G0 with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, u, t_cut
        real(kind=8) :: dME_G0_REG


        ! dME_G0_REG = 9D0/2D0 * ( 3D0 - u*t/s**2 &
        dME_G0_REG = 9D0/2D0 * ( 3D0 - u*t/(s + t_cut)**2 &
                   - s*u/(t - t_cut)**2 - t*s/(u - t_cut)**2 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G1( s, t, u, m1, m2 )
!!      The Matrix Element component G1.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m1, m2
        real(kind=8) :: dME_G1


        dME_G1 = 2D0/9D0 * 2D0/t**2 * ( (s - (m1**2 + m2**2) )**2 &
               + (m1**2 + m2**2 - u)**2 + 2D0*t*(m1**2 + m2**2) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G1_REG( s, t, u, m1, m2, t_cut )
!!      The Matrix Element component G1(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, u, m1, m2, t_cut
        real(kind=8) :: dME_G1_REG


        dME_G1_REG = 2D0/9D0 * 2D0/(t - t_cut)**2    &
                   * ( (s - (m1**2 + m2**2) )**2     &
                   + (m1**2 + m2**2 - u)**2 + 2D0*t*(m1**2 + m2**2) )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G2( s, t, u, m )
!!      The Matrix Element component G2.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m
        real(kind=8) :: dME_G2


        dME_G2 = - 2D0/27D0 * 4D0/t/u * (s - 2D0*m**2) * (s - 6D0*m**2)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G2_REG( s, t, u, m, t_cut, u_cut )
!!      The Matrix Element component G2.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m, t_cut, u_cut
        real(kind=8) :: dME_G2_REG


        dME_G2_REG = - 2D0/27D0 * 4D0 / (t - t_cut) / (u - u_cut) &
                   * (s - 2D0*m**2) * (s - 6D0*m**2)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G3( s, t, u, m )
!!      The Matrix Element component G3.
!       Cf. V. Borchers et al., Phys.Rev.C 62 (2000) 064903.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m
        real(kind=8) :: dME_G3


        dME_G3 = 3D0/16D0 * 4D0 * (m**2 - t) * (m**2 - u) / s**2          &
               + 1D0/12D0 * ( 2D0 * (m**2 - u) * (m**2 - t)               &
                - 4D0 * m**2 * (m**2 + t) ) / (t - m**2)**2               &
               + 1D0/12D0 * ( 2D0 * (m**2 - t) * (m**2 - u)               &
                - 4D0 * m**2 * (m**2 + u) ) / (u - m**2)**2               &
               - 1D0/96D0 * 4D0 * m**2 * (s - 4D0*m**2)                   &
                / (m**2 - t) / (m**2 - u)                                 &
               + 3D0/32D0 * 4D0 * ( t * (s + t) - m**4 ) / s / (m**2 - t) &
               + 3D0/32D0 * 4D0 * ( u * (s + u) - m**4 ) / s / (m**2 - u)


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_G3_REG( s, t, u, m, s_cut )
!!      The Matrix Element component G3.
        implicit none
        real(kind=8), intent(in) :: s, t, u, m, s_cut
        real(kind=8) :: dME_G3_REG


        dME_G3_REG = 3D0/16D0 * 4D0 * (m**2 - t) * (m**2 - u) / (s + s_cut)**2 &
                   + 1D0/12D0 * ( 2D0 * (m**2 - u) * (m**2 - t)                &
                    - 4D0 * m**2 * (m**2 + t) ) / (t - m**2)**2                &
                   + 1D0/12D0 * ( 2D0 * (m**2 - t) * (m**2 - u)                &
                    - 4D0 * m**2 * (m**2 + u) ) / (u - m**2)**2                &
                   - 1D0/96D0 * 4D0 * m**2 * (s - 4D0*m**2)                    &
                     / (m**2 - t) / (m**2 - u)                                 &
                   + 3D0/32D0 * 4D0 * ( t * (s + t) - m**4 ) &
                    / ( s + s_cut ) / ( m**2 - t ) &
                   + 3D0/32D0 * 4D0 * ( u * (s + u) - m**4 ) &
                    / ( s + s_cut ) / ( m**2 - u )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF0( s, t )
!!      The primitive function of the Matrix Element component G0.
        implicit none
        real(kind=8), intent(in) :: s, t
        real(kind=8) :: dME_PF0


        dME_PF0 = 3D0*t + t**2/2D0/s + t**3/3D0/s**2 - s**2/t &
                + s*LOG( ABS( t / (s + t) ) ) - s**2/(s + t)
        dME_PF0 = 9D0 / 2D0 * dME_PF0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF0_REG( s, t, t_cut )
!!      The primitive function of the Matrix Element component G0 w/ regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, t_cut
        real(kind=8) :: dME_PF0_REG


        dME_PF0_REG = 3D0*t + s*t**2/2D0/(s + t_cut)**2             &
                    + t**3/3D0/(s + t_cut)**2                       &
                    - s**2/(t - t_cut)                              &
                    + s*LOG( ABS( (t - t_cut) / (s + t + t_cut) ) ) &
                    - s*t_cut/(t - t_cut) - (s**2 + s*t_cut)/(s + t + t_cut)
        dME_PF0_REG = 9D0 / 2D0 * dME_PF0_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_stu( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_stu, M


        M = m1**2 + m2**2
        dME_PF1_stu = - 2D0*(s - M)**2 / t &
                    + 2D0 * S * LOG( ABS(t) ) + t
        dME_PF1_stu = 4D0 / 9D0 * dME_PF1_stu


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_stu_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_stu_REG, M


        M = m1**2 + m2**2
        dME_PF1_stu_REG = - 2D0*(s - M)**2 / (t - t_cut)  &
                        - 2D0*(s - M)*t_cut / (t - t_cut) &
                        - t_cut**2 / (t - t_cut)          &
                        + (t - t_cut)                     &
                        + (2D0*s + 2D0*t_cut) * LOG( ABS(t - t_cut) )
        dME_PF1_stu_REG = 4D0 / 9D0 * dME_PF1_stu_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_sut( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(s,u,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_sut, M


        M = m1**2 + m2**2
        dME_PF1_sut = 2D0*(s - M)**2 / (2D0*M - s - t) &
                    - (2D0*M - s - t)                  &
                    - 2D0*s*LOG( ABS(2D0*M - s - t) )
        dME_PF1_sut = 4D0 / 9D0 * dME_PF1_sut


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_sut_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,u,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_sut_REG, M


        M = m1**2 + m2**2
        dME_PF1_sut_REG = ( (s - M)**2 + (s - M + t_cut)**2 ) &
                          / (2D0*M - s - t - t_cut)           &
                        - (2D0*M - s - t - t_cut)             &
                        - 2D0*(s + t_cut)*LOG( ABS(2D0*M - s - t - t_cut) )
        dME_PF1_sut_REG = 4D0 / 9D0 * dME_PF1_sut_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_ust( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(u,s,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_ust, M


        M = m1**2 + m2**2
        dME_PF1_ust = 2D0 * (t - M)**3 / 3D0 / s**2 &
                    + t**2 / s + t
        dME_PF1_ust = 4D0 / 9D0 * dME_PF1_ust


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_ust_REG( s, t, m1, m2, s_cut )
!!      The primitive function of the Matrix Element component
!!       G1(s,u,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, s_cut
        real(kind=8) :: dME_PF1_ust_REG, M


        M = m1**2 + m2**2
        dME_PF1_ust_REG = 2D0 * (t - M)**3 / 3D0 / (s + s_cut)**2 &
                        + t**2 / (s + s_cut) + s**2 * t / (s + s_cut)**2
        dME_PF1_ust_REG = 4D0 / 9D0 * dME_PF1_ust_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_uts( s, t, m1, m2 )
!!      The primitive function of the Matrix Element component G1(u,t,s).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2
        real(kind=8) :: dME_PF1_uts, M


        M = m1**2 + m2**2
        dME_PF1_uts = -2D0*(M - s)**2 / t + 2D0*s*LOG( ABS(t) ) + t
        dME_PF1_uts = 4D0 / 9D0 * dME_PF1_uts


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF1_uts_REG( s, t, m1, m2, t_cut )
!!      The primitive function of the Matrix Element component
!!       G1(u,t,s) with the regulator
        implicit none
        real(kind=8), intent(in) :: s, t, m1, m2, t_cut
        real(kind=8) :: dME_PF1_uts_REG, M


        M = m1**2 + m2**2
        dME_PF1_uts_REG = -( (M - s)**2 + (M - s - t_cut)**2 ) / (t - t_cut) &
                        + 2D0*(s + t_cut)*LOG( ABS(t - t_cut) )              &
                        + (t - t_cut)
        dME_PF1_uts_REG = 4D0 / 9D0 * dME_PF1_uts_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_stu( s, t, m1 )
!!      The primitive function of the Matrix Element component G2(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF2_stu, M


        M = m1**2 + m1**2
        dME_PF2_stu = (s - M) * (s - 3D0*M) &
                    * LOG( ABS( t / (2D0*M - s - t) ) ) / (2D0*M - s)
        dME_PF2_stu = - 8D0 / 27D0 * dME_PF2_stu


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_stu_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component
!!       G2(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF2_stu_REG, M


        M = m1**2 + m1**2
        dME_PF2_stu_REG = (s - M) * (s - 3D0*M) &
                        * LOG( ABS( (t - t_cut) / (2D0*M - s - t - t_cut) ) ) &
                        / (2D0*M - s - 2D0*t_cut)
        dME_PF2_stu_REG = - 8D0 / 27D0 * dME_PF2_stu_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_ust( s, t, m1 )
!!      The primitive function of the Matrix Element component G2(u,s,t).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF2_ust, M


        M = m1**2 + m1**2
        dME_PF2_ust = 2D0*t + t**2/2D0/s &
                    + (s - M**2/s) * LOG( ABS(t) )
        dME_PF2_ust = - 8D0 / 27D0 * dME_PF2_ust


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF2_ust_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component
!!       G2(u,s,t) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF2_ust_REG, M


        M = m1**2 + m1**2
        dME_PF2_ust_REG = ( 2D0*(s + t_cut)*(t - t_cut)          &
                        + 1D0/2D0*(t - t_cut)**2                 &
                        + (s**2 - M**2 + 2D0*s*t_cut + t_cut**2) &
                        * LOG( ABS(t - t_cut) ) )                &
                        / (s + t_cut)
        dME_PF2_ust_REG = - 8D0 / 27D0 * dME_PF2_ust_REG


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_stu( s, t, m1 )
!!      The primitive function of the Matrix Element component G3(s,t,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF3_stu, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * ( M*t**2 - s*t**2/2D0 &
           - t**3/3D0 - M**2*t + M*s*t ) / s**2
        T2 = 2D0 * ( - t - s*LOG( ABS(M - t) ) ) &
           - 4D0 * M * ( 2D0*M/(M - t) + LOG( ABS(M - t) ) )
        T3 = 2D0 * ( s*LOG( ABS(s + t - M) ) - t ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M*(s - 4D0*M)/s * LOG( ABS( (s + t - M) / (M - t) ) )
        T5 = 4D0 * ( s*M - s*t + 2D0*M**2 - 2D0*M*t &
           - (M - t)**2 / 2D0 - s*M*LOG( ABS(M - t) ) ) / s
        T6 = 4D0 * ( (s + t - M)**2 / 2D0 - s**2 - s*t - s*M - 2D0*M*t &
           + 2D0*M**2 + s*M*LOG( ABS(s + t - M) ) ) / s
        ! T1 = 0D0
        dME_PF3_stu = 3D0/16D0*T1 + 1D0/12D0*T2 &
                    + 1D0/12D0*T3 - 1D0/96D0*T4 &
                    + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_stu_REG( s, t, m1, s_cut )
!!      The primitive function of the Matrix Element component
!!       G3(s,t,u) with the regulator.
        implicit none
        real(kind=8), intent(in) :: s, t, m1, s_cut
        real(kind=8) :: dME_PF3_stu_REG, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * ( M*t**2 - s*t**2/2D0 &
           - t**3/3D0 - M**2*t + M*s*t ) / (s + s_cut)**2
        T2 = 2D0 * ( - t - s*LOG( ABS(M - t) ) ) &
           - 4D0 * M * ( 2D0*M/(M - t) + LOG( ABS(M - t) ) )
        T3 = 2D0 * ( s*LOG( ABS(s + t - M) ) - t ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M*(s - 4D0*M)/s * LOG( ABS( (s + t - M) / (M - t) ) )
        T5 = 4D0 * ( s*M - s*t + 2D0*M**2 - 2D0*M*t &
           - (M - t)**2 / 2D0 - s*M*LOG( ABS(M - t) ) ) / (s + s_cut)
        T6 = 4D0 * ( (s + t - M)**2 / 2D0 - s**2 - s*t - s*M - 2D0*M*t &
           + 2D0*M**2 + s*M*LOG( ABS(s + t - M) ) ) / (s + s_cut)
        dME_PF3_stu_REG = 3D0/16D0*T1 + 1D0/12D0*T2 &
                        + 1D0/12D0*T3 - 1D0/96D0*T4 &
                        + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_tsu( s, t, m1 )
!!      The primitive function of the Matrix Element component G3(t,s,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1
        real(kind=8) :: dME_PF3_tsu, T1, T2, T3, T4, T5, T6, M


        M  = m1**2
        T1 = 4D0 * (M - s) * ( M/t - s/t + LOG( ABS(t) ) )
        T2 = 2D0 * ( s*t + t**2/2D0 - M*t ) / (M - s) &
           - 4D0 * M * (M + s) / (s - M)**2 * t
        T3 = 2D0 * (M - s) * LOG( ABS(s + t - M) ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M / (M - s) * ( s + t - M &
           - s*LOG( ABS(s + t - M) ) - 3D0*M*LOG( ABS(s + t - M) ) )
        T5 = 4D0 * ( s*t + (s**2 - M**2) * LOG( ABS(t) ) ) / (M - s)
        T6 = 4D0 * ( (3D0*M**2 - 4D0*M*s + s**2) / (s - M) &
           * LOG( ABS( t / ( s + t - M ) ) ) &
           - (2D0*M - s) * LOG( ABS(s + t - M) ) )
        dME_PF3_tsu = 3D0/16D0*T1 + 1D0/12D0*T2 &
                    + 1D0/12D0*T3 - 1D0/96D0*T4 &
                    + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        function dME_PF3_tsu_REG( s, t, m1, t_cut )
!!      The primitive function of the Matrix Element component G3(t,s,u).
        implicit none
        real(kind=8), intent(in) :: s, t, m1, t_cut
        real(kind=8) :: dME_PF3_tsu_REG, T1, T2, T3, T4, T5, T6, A, B, C, D, M


        M  = m1**2
        T1 = 4D0 * (M - s) * ( s/(t_cut - t) + t_cut/(t_cut - t) &
           + LOG( ABS(t - t_cut) ) + M/(t - t_cut) )
        T2 = 2D0 * ( s*t + t**2/2D0 - M*t ) / (M - s) &
           - 4D0 * M * (M + s) / (s - M)**2 * t
        T3 = 2D0 * (M - s) * LOG( ABS(s + t - M) ) &
           - 4D0 * M * ( 2D0*M / (M - s - t) - LOG( ABS(M - s - t) ) )
        T4 = 4D0 * M / (M - s) * ( s + t - M &
           - s*LOG( ABS(s + t - M) ) - 3D0*M*LOG( ABS(s + t - M) ) )
        T5 = 4D0 * ( -M**2*LOG( ABS(t - t_cut) ) &
           + s * ( (s + t_cut)*LOG( ABS(t - t_cut) ) + (t - t_cut) ) ) &
           / (M - s)
         A = ( 3D0*M**2 - 4D0*M*s + s**2 ) / ( s + t_cut - M )
         B = - A
         C = - t_cut * ( 2D0*M - s ) / ( s + t_cut - M )
         D = - ( 2D0*M - s ) - C
        T6 = 4D0 * ( (A + C) * LOG( ABS(t - t_cut) ) &
           + (B + D) * LOG( ABS(s + t - M) ) )
        dME_PF3_tsu_REG = 3D0/16D0*T1 + 1D0/12D0*T2 &
                        + 1D0/12D0*T3 - 1D0/96D0*T4 &
                        + 3D0/32D0*T5 + 3D0/32D0*T6


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine find_par( icp, tcp )
!!      Finds out the binary collision (icp) with minimum colli. time (tcp).
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MCLIS=280000)
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol


        icp = MINLOC( tc( 1, 1:icol ), DIM=1 )
        tcp = MINVAL( tc( 1, 1:icol ), DIM=1 )


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine collis( icp )
!!      Performs parton-parton collision & updates particle list
!       icp: the current icp-th collision pair
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        COMMON/PYJET2/I_PTR2AA(KSZJ),K6K7K8(KSZJ,3),P6P7(KSZJ,2)
        COMMON/PYJET3/ N00_PY8, N00_DIQ, I_COL_MAX(2)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/shower_1/ KSH(14,10), PSH(14,10), VSH(14,10)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_deadcone,i_LPM,i_diquark, &
                    ipad25,para1_1,para1_2
        common/syspar_p/rsig1,pio,tcut
        common/ctllist_p/nreac(100),nrel
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        common/work7/reac(100),crose(100)
!       lc(1,i), lc(2,i): line numbers in the particle list of i-th coll. pair
!       lc(3,i), lc(4,i): flavor codes of scattered particles of i-th coll. pair
!       lc(5,i): the internal process number
!       lc(6,i): event type. =1, elastic scattering; =2, inelastic; =3, shower.
!       tc: the collision time of the colliding pair
!       sig_tsmp_Cthetas_p(1,i): the total cross section of the colliding pair
!       sig_tsmp_Cthetas_p(2,i): the sampled mom. transfer "tsmp" of the pair
!       sig_tsmp_Cthetas_p(3,i): the scattering angle COS(theta_s) of the pair
!       sig_tsmp_Cthetas_p(4,i): the momentum modulus of th pair in the CM frame
!       ic, jc: line number of colliding particles
!       kf3, kf4: kf codes of the collided pair after the collision
!       nreac(i): statistics of the # of successful i-th collision
!       nrel: statistics of the # of collision blocked
!       nsca: number of particles after collision (+shower)
        integer, intent(in) :: icp
        dimension pi(4), pj(4)
        dimension pij(4), b(3)
        real(kind=8) :: Lambda_QCD2


!       Gets contents of the current collision pair.
        ic     = lc(1,icp)
        jc     = lc(2,icp)
        ic_sh  = 0
        jc_sh  = 0
        kf3    = lc(3,icp)
        kf4    = lc(4,icp)
        i_proc = lc(5,icp)
        sig         = sig_tsmp_Cthetas_p(1,icp)
        tsmp        = sig_tsmp_Cthetas_p(2,icp)
        cos_theta_s = sig_tsmp_Cthetas_p(3,icp)
        p_modulus   = sig_tsmp_Cthetas_p(4,icp)

!       kf1 and kf2: kf code of the colliding pair before interaction
        kf1 = K(ic,2)
        kf2 = K(jc,2)
!       Selection of the form of differential cross section in the LO-pQCD.
        i_LOpQCD = INT( adj1(20) )


!-------------------------------------------------------------------------------
!----------------------------   Paticle Collision   ----------------------------
!       Number of scattered (+ showered) particles.
        nsca  = 0
        pi(4) = P(ic,4)
        pj(4) = P(jc,4)
        do i=1,3,1
            pi(i) = P(ic,i)
            pj(i) = P(jc,i)
            b(i)  = ( pi(i) + pj(i) ) / ( pi(4) + pj(4) )
        end do
        pij = pi + pj
!       Invariant mass squared of the colliding pair.
        eiej2 = pij(4)*pij(4) - pij(1)*pij(1) - pij(2)*pij(2) - pij(3)*pij(3)
!       Boosts to the current CMS of the colliding pair.
        call lorntz( 0, b, pi, pj )
        dm1 = P(ic,5)
        dm2 = P(jc,5)
        dm3 = dm1
        dm4 = dm2
        am1 = dm1
        am2 = dm2
        am3 = dm3
        am4 = dm4
        i_proc1 = MOD( i_proc, 10 )
        ! Mass treatment corresponds to "fsig".
        select case( i_proc1 )
        ! Inelstic
        case( 2, 6, 7 )
            dm3 = amass( kf3 )
            dm4 = amass( kf4 )
            am3 = dm3
            am4 = dm4
        ! Elastic
        case default
!       Small angle and zero-mass of light quarks approximations for el. scat.
            if( i_LOpQCD == 1 .OR. i_LOpQCD == 3 )then
                if( ABS(kf1) < 4 ) am1 = 0D0
                if( ABS(kf2) < 4 ) am2 = 0D0
                am3 = am1
                am4 = am2
            end if
        end select

!Lei_debug
!       Special sampling sunroutine for massive processes. Function 'samplr_t'
!        just gives the simple sapmling without masses.
        if( am1+am2+am3+am4 > 1D-5 )then
            t_min = sig_tsmp_Cthetas_p(2,icp)
            t_max = sig_tsmp_Cthetas_p(3,icp)
            call sample_t_massive( t_min, t_max, am1, am2, am3, am4, eiej2, &
                                   i_proc1, i_LOpQCD, icp )
            tsmp        = sig_tsmp_Cthetas_p(2,icp)
            cos_theta_s = sig_tsmp_Cthetas_p(3,icp)
        end if
!Lei_debug

!       Calculates the direction cosines of momentum of one colliding
!        particle after scattering relative to the momentum
!        of the corresponding particle before scattering.
        Ctheta_s = cos_theta_s
!       Random 2*pi azimuthal angle.
        phi_s = 2D0 * pio * PYR(1)
!       Rescales if Ctheta_s overflows.
        if( ABS(Ctheta_s) > 1D0 )then
            Ctheta_s = Ctheta_s / ABS( Ctheta_s ) - 1D-10
        end if
        Stheta_s = SQRT( 1D0 - Ctheta_s*Ctheta_s )
        Cphi_s   = COS( phi_s )
        Sphi_s   = SIN( phi_s )

!       pi and pj, as input into "rotate", are four momentum of the colliding
!        pair (before scattering and rotation), as output are four
!        momentum of that pair after scattering and rotation
!       p_modulus: momentum modulus of pi or pj, both are equal in their CMS.
        call rotate( Ctheta_s, Stheta_s, Cphi_s, Sphi_s, p_modulus, pi, pj )
        pi(4) = SQRT( pi(1)*pi(1) + pi(2)*pi(2) + pi(3)*pi(3) + dm3*dm3 )
        pj(4) = SQRT( pj(1)*pj(1) + pj(2)*pj(2) + pj(3)*pj(3) + dm4*dm4 )

!       Boosts back to Lab.
        call lorntz( 1, b, pi, pj )

!       Updates particle list after the scattering (note: line numbers of the
!        colliding pair after scattering are the same as ones before scattering)
        ! ic
        K(ic,2) = kf3
        do i=1,4,1
            P(ic,i) = pi(i)
        end do
!       For the inelastic scattering.
        if( kf3 /= kf1 ) P(ic,5) = dm3
        ! jc
        K(jc,2) = kf4
        do i=1,4,1
            P(jc,i) = pj(i)
        end do
        if( kf4 /= kf2 ) P(jc,5) = dm4
        nsca = 2

!       Counts cross sections.
        reac(i_proc)  = reac(i_proc)  + 1D0
        crose(i_proc) = crose(i_proc) + sig
        nreac(i_proc) = nreac(i_proc) + 1
!----------------------------   Paticle Collision   ----------------------------
!-------------------------------------------------------------------------------


!       Without the time-like shower.
        KSH = 0
        PSH = 0D0
        VSH = 0D0
        if( i_time_shower == 0 ) return


!-------------------------------------------------------------------------------
!-------------------------   Time-like Shower Finding  -------------------------
!       Mandelstam variables.
        xs  = eiej2
        xt  = tsmp
        xu  = am1*am1 + am2*am2 + am3*am3 + am4*am4 - xs - xt
!       Squareed momentum transfer (Q^2) to be used in the time-like shower.
        ! pT^2, same as that in PYTHIA ( Q^2_{hard} ).
        ! Massless.
        ! qc2 = xt * xu / xs
        ! Massive.
        qc2 = ( xt * xu - am1*am1*am2*am2 - am3*am3*am4*am4 ) / xs
        ! Corresponding to 4*tl0.
        ! qc2 = 4D0 * qc2
        ! qc2 = MIN( -xt, -xu )
!       Effective coupling constant for shower.
        ! Constant alpha_s.
        as  = adj1(2)
        as0 = adj1(2)
        ! Renormalization scale mu_R = Q.
        ! Q2 = s - dm1*dm1 - dm2*dm2
        Q2 = qc2
        ! Number of active quarks.
        Nf = 3
        if( Q2 > amass(4)*amass(4) ) Nf = 4
        if( Q2 > amass(5)*amass(5) ) Nf = 5
        ! QCD scale.
        Lambda_QCD2 = adj1(25)*adj1(25)
        ! Running alpha_s.
        if( as0 < 1D-10 ) as = func_alpha_s( Nf, Q2, Lambda_QCD2 )
        if( as  < 1D-10 ) return

!       Prepares for the time-like shower.

        ! Final state of the scattering (nsca=2).
        ! Also the initial state of the time-like branching.
        !
        ! This part will be used for the caluculation of the shower subroutine.
        !
        ! pip(1-3,*): momentum of the particle after collision/shower
        ! pip(4,*): energy
        ! pip(5,*): virtuality
        ! pip(6,*): x value (energy fraction)
        ! kpip(*): flavor code
        kpip(1) = kf3
        kpip(2) = kf4
        do i=1,4,1
            pip(i,1) = pi(i)
            pip(i,2) = pj(i)
        end do
        pip(5,1) = qc2
        pip(5,2) = qc2
        pip(6,1) = 1D0
        pip(6,2) = 1D0
        !
        ! This part will be used to update the particle and collision time list.
        !
        !  1-2 : the initial state of the shower (final state of scattering)
        !  3-5 : two radiators with shifted momentum or two 1st shower partons
        !  4-6 : two showerd partons of '3' and '5'
        !
        !  KSH1: status, K2: id, K3: mother1, K4: daughter1, K5: daughter2,
        !  KSH6: mother2, K7: color tag, K8: anti-color tag, K9: shower type.
        !  KSH9: =31, q -> qg; =32, g -> gg; =33, g -> qqbar; (light)
        !        =41, c -> cg;             ; =43, g -> ccbar; (heavy c)
        !        =51, b -> bg;             ; =53, g -> bbbar. (heavy b)
        !  Shower case.
        !        = 1, 1+2 radiation, two recoilers.
        !        = 2, parton 1 radiates, recoils paron 2.
        !        = 3, parton 2 radiates, recoils paron 1.
        !
        !  PSH1-PSH5: px, py, pz, e, m; PSH6: scale(Q2), PSH7: spin/polarization
        !  VSH1-VSH5: x, y, z, t, tau(formation time in the Lab frame)
        ! 1
        KSH(1,1) = K(ic,1)
        KSH(1,2) = K(ic,2)
        KSH(1,3) = ic - 1
        KSH(1,7) = K6K7K8(ic,2)
        KSH(1,8) = K6K7K8(ic,3)
        do i=1,5,1
            PSH(1,i) = P(ic,i)
            VSH(1,i) = V(ic,i)
        end do
        PSH(1,6) = P6P7(ic,1)
        PSH(1,7) = P6P7(ic,2)
        ! 2
        KSH(2,1) = K(jc,1)
        KSH(2,2) = K(jc,2)
        KSH(2,3) = jc - 1
        KSH(2,7) = K6K7K8(jc,2)
        KSH(2,8) = K6K7K8(jc,3)
        do i=1,5,1
            PSH(2,i) = P(jc,i)
            VSH(2,i) = V(jc,i)
        end do
        PSH(2,6) = P6P7(jc,1)
        PSH(2,7) = P6P7(jc,2)

!       n_call_tili: statistics of the times calling "ti_li1".
        n_call_tili = 0
!       tl0: the virtuality cutoff of the time-like branching, i.e. Mu0^2 (Q0^2).
        tl0 = adj1(24)

!       The time-like branching trial. "nsca" may increase in "ti_li1".
        nsca0 = nsca
        do i = 1, nsca0, 1
            kfk = kpip(i)
            if( ABS(kfk) > 21 ) cycle
            ! Phase space for branching vanishes at a virtuality of 4*tl0.
            !  cf. Eq.25 in B.R. Webber, Ann. Rev. Nucl. Part. Sci. 36 (86) 253.
            if( pip(5,i) > 4D0*tl0 )then
                ea = pip(4,i)
                qa = pip(5,i)
                do j=1,3
                    pa(j) = pip(j,i)
                end do
                xa = pip(6,i)
                ! Core subroutine. Performs the branching.
                call ti_li1( n_call_tili, i, i_succ )
                n_call_tili = n_call_tili + 1
                ! No showers occurred. Takes a carbon copy for the radiator.
                if( i_succ == 0 )then
                    do j=1,10,1
                        KSH( 2*i+1, j ) = KSH(i,j)
                        PSH( 2*i+1, j ) = PSH(i,j)
                        VSH( 2*i+1, j ) = VSH(i,j)
                    end do
                ! Labels the shower type to radiator.
                else
                    KSH(i,9) = KSH( 2*i+1, 9 )
                end if
            end if
        end do
!-------------------------   Time-like Shower Finding  -------------------------
!-------------------------------------------------------------------------------


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine sample_t_massive( t_min, t_max, m1, m2, m3, m4, s, &
                                     i_proc, i_LOpQCD, icp )
!!      Special sampling sunroutine for massive processes.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MCLIS=280000)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        integer, intent(in) :: i_proc, i_LOpQCD, icp
        real(kind=8), intent(in) :: t_min, t_max, m1, m2, m3, m4
        real(kind=8) :: Lambda_QCD2


!       Resets.
        t = 0D0
        cos_theta_s = 0D0
        sig_tsmp_Cthetas_p(2,icp) = 0D0
        sig_tsmp_Cthetas_p(3,icp) = 0D0

!       Effective coupling constant.
        ! Constant alpha_s.
        alpha_s  = adj1(2)
        alpha_s0 = adj1(2)
        ! Renormalization scale mu_R = Q.
        Q2 = s - m1*m1 - m2*m2
        ! Number of active quarks.
        Nf = 3
        if( Q2 > amass(4)*amass(4) ) Nf = 4
        if( Q2 > amass(5)*amass(5) ) Nf = 5
        ! QCD scale.
        Lambda_QCD2 = adj1(25)*adj1(25)
        ! Running alpha_s.
        if( alpha_s0 < 1D-10 ) alpha_s = func_alpha_s( Nf, Q2, Lambda_QCD2 )
        if( alpha_s  < 1D-10 ) return
!       t_cut: the cutoff used to regulate the 't'/'u' divergence.
        ! Constant cutoff (screening masses).
        t_cut  = adj1(3)
        t_cut0 = adj1(3)
        ! Gluon screening mass.
        dmD2   = t_cut
        ! Quark medium mass.
        dmq2   = t_cut
        ! Dynamic screening masses.
        if( t_cut0 < 1D-10 )then
            T_effictive = ABS( adj1(3) )
            ! Boltzmann or Bose-Einstein & Fermi-Dirac statistics in screening.
            ! i_statistic = 0
            i_statistic = 1
            ! Gluon screening mass.
            dmD2  = func_mD2( alpha_s, Nf, T_effictive, i_statistic )
            t_cut = dmD2
            ! Quark medium mass.
            dmq2  = func_mq2( alpha_s, T_effictive, i_statistic )
        end if

!       Full LO-pQCD forms for el. and inel. cross sections. Rejection sampling.
        do while(.true.)
            rand_num_1 = PYR(1)
            t = t_min + rand_num_1 * ( t_max - t_min )
            dsigma_dt     = dsigma_dt_func( s, t, t_cut, m1, m2,    &
                                            m3, m4, i_proc, i_LOpQCD )
            dsigma_dt_max = MAX(                                    &
                        dsigma_dt_func( s, t_min, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ), &
                        dsigma_dt_func( s, t_max, t_cut, m1, m2,    &
                                        m3, m4, i_proc, i_LOpQCD ) )
            prob_accept = dsigma_dt / dsigma_dt_max
            rand_num_2  = PYR(1)
            if( rand_num_2 <= prob_accept ) exit
        end do

!       The cosine of the scattering angle theta_s.
        ! cos_theta_s = 1D0 - 2D0 * ( t_0 - tsmp ) / ( t_0 - t_pi )
!       Equivalently.
        cos_theta_s = 1D0 + 2D0*s*t / ( s - m1 - m2 ) / ( s - m1 + m2 )

!       Feeds back.
        sig_tsmp_Cthetas_p(2,icp) = t
        sig_tsmp_Cthetas_p(3,icp) = cos_theta_s


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine ti_li1( n_call_tili, ij, i_succ )
!!      Performs the time-like branching (along main chain until
!!        tl0=Mu0^2 (Q0^2) for the branching parton).
!       n_call_tili: statistics the times calling "ti_li1"
!       ij: order # of spliting parton
!       i_succ: successful flag
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (pio=3.141592653589793D0, hbarc=0.197327D0)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/shower_1/ KSH(14,10), PSH(14,10), VSH(14,10)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/sa25/i_inel_proc,i_time_shower,i_deadcone,i_LPM,i_diquark, &
                    ipad25,para1_1,para1_2
        dimension pa1(4),pa2(4),p00(4)
        integer :: n_sample_q


!       i_stop=0 means the branching will go on, i_stop=1 means "stop".
        i_stop = 0
        i_succ = 0
        pa1 = 0D0
        pa2 = 0D0
        p00 = 0D0

        n_sample_q = 0
10      continue
        n_sample_q = n_sample_q + 1

!       Virtualities in time-like branchs are all the mass squared (m**2).
!        and m**2 = E**2 - |p|**2 > 0
!------------------------------------------------------------------
!       Branching:   A  -> A1   +   A2
!       Vituality:  qqb    qqa     qqc
!       ID No.   :  kfk    kf1     kf2
!        x       :  xa x'  xa1 x   xa2
!       Momentum :  pa     pa1     pa2
!
!       Equation : qqb = qqa/z + qqc/(1-z) + pT2 / ( z*(1 - z) ) (*)
!-------------------------------------------------------------------
!       Above Eq. is deduced from the relation of mass and energy (Einstein
!        relation) and four-momentum conservation, cf. note made in
!        Webber's paper.
!-------------------------------------------------------------------
!       A2 will be always the radiated one.

!-------Flavors of A1, A2 will be decided in the ensuing module.
        q2max = qa
        kf0   = kfk
        call suda_ti( i_fail, kf0, i_split0, q2max, q2, z )
!       Time-like Sudakov factor, input: kf0 & q2max, output: q2 & z.
!       i_fail: =1, failed in sampling q2 and z
!               =0, succeeded in sampling q2 and z

!       Forcibly stops.
        if( i_fail == 1 ) then
            qa = tl0
            return
        end if

        qqb = q2
        if( n_call_tili > 2 ) qqb = qa
        zab = z
!       "zab" is the sampled ratio of x/x' in this branching.
        if( ABS(zab) <= 1D-10 ) return
        xa1 = xa*zab
        xa2 = xa - xa1
        kf1 = 0
        kf2 = 0
!       Initial state is quark.
        if( kfk /= 21 ) then
        ! Branching q -> qg.
            kf1 = kfk
            kf2 = 21
!       Initial state is gluon.
        else
            ! Branching g -> gg.
            if( i_split0 == 2 ) then
                kf1 = 21
                kf2 = 21
            ! Branching g -> qqbar
            else
                ! Only consider g -> gg now.
                IF( .FALSE. )THEN
                    ! ea = pip(4,ij): energy of spliting particle ij
                    ! pip(1-3,ij): three momentum of spliting particle ij
                    eaa = ea*ea - pip(1,ij)*pip(1,ij) - pip(2,ij)*pip(2,ij) &
                        - pip(3,ij)*pip(3,ij)
                    if( eaa >= 4D0*amass(2)*amass(2) )then
                        ! Invariant mass of the initial state gluon.
                        eaa = SQRT(eaa)
                        call break_f( eaa, kf_out, dm_out )
                        kf1 =  kf_out
                        kf2 = -kf1
                    end if
                END IF
            end if
        end if
!----------Finished.

        n_sample_pT = 0
300     continue
        n_sample_pT = n_sample_pT + 1
        if( n_sample_pT > 1000 ) goto 10

        q2max = zab*( qqb - tl0 / (1D0 - zab) )
!       The max. value of A1's virtuality will be decided by the eqution(*)
!         with qqc=tl0 & pT2=0 .
        call suda_ti( i_fail, kf1, i_split, q2max, q2, z )
        if( i_fail == 1 ) then
!       qqa=tl0 means the forward evolution has finished, thus set i_stop=1.
            qqa = tl0
            i_stop = 1
        else
            qqa = q2
        end if
        q2max = ( qqb - qqa/zab ) * ( 1D0 - zab )
!       The max. value of A2's virtuality will be decided by the Eq(*) w/ pT2=0.
        call suda_ti( i_fail, kf2, i_split, q2max, q2, z )
        if( i_fail == 1 ) then
            qqc = tl0
        else
            qqc = q2
        endif
!       pT2: the squared transverse momentum determined by the Eq.(*)
        pT2 = ( qqb - qqa/zab - qqc/(1D0 - zab) ) * ( 1D0 - zab ) * zab

!0130   The following block gives the momentum of A1, A2 i.e. pa1, pa2
!        according to the definition of z & pT.
        pa02 = pa(1)*pa(1) + pa(2)*pa(2) + pa(3)*pa(3)
        ! Rotates to the frame in which the momentum of A is along +Z direction.
        ! Third momentum of A. pT = 0.
        paz  = SQRT(pa02)
        ! Third momentum of A1.
        pz1  = zab*paz
        ! Third momentum of A2.
        pz2  = ( 1D0 - zab ) * paz
        ! Momentum modulus of A1.
        ppa1 = SQRT( pz1*pz1 + pT2 )
        ! Momentum modulus of A2.
        ppa2 = SQRT( pz2*pz2 + pT2 )
        ppt  = SQRT(pT2)

!       Imposes the dead-cone effect.
        theta_split = ATAN( ppt / pz2 )
        E_mother    = PSH(ij,4)
        dm_mother   = PSH(ij,5)
        theta_dead  = dm_mother / E_mother

!       Suppresses the mass-related small angle emmision. Resamples.
        if( n_sample_q > 10 ) return
        if( i_deadcone /= 0 .AND. theta_split < theta_dead ) goto 300

!       Calculates the formation time of a radiated gluon (quark?)
!        in the Lab frame of the scatted (radiators) partons for the LPM effect.
!        t_form = 2*E / ( pT^2 + x^2 * M^2 ), M: mass of the mother (radiator)
        if( i_LPM == 0 )then
            t_form1 = 1D-8
            t_form2 = 1D-8
        else
            dm1 = amass(kf1)
            if( kf1 == KSH(ij,2) ) dm1 = dm_mother
            dm2 = amass(kf2)
            if( kf2 == KSH(ij,2) ) dm2 = dm_mother
            E_split1 = SQRT( pz1*pz1 + pT2 + dm1*dm1 )
            t_form1  = hbarc * 2D0*E_split1 / ( pT2 &
                     + zab * zab * dm_mother * dm_mother )
            E_split2 = SQRT( pz2*pz2 + pT2 + dm2*dm2 )
            t_form2  = hbarc * 2D0*E_split2 / ( pT2  &
                     + ( 1D0 - zab ) * ( 1D0 - zab ) * dm_mother * dm_mother )
        end if

        sctas1 = ppt / ppa1
        cctas1 = pz1 / ppa1
        phi = 2D0*pio*PYR(1)
        cfis = COS(phi)
        sfis = SIN(phi)
!       Direction cosines of A1 relative to momentum of A
!       Note : transverse momentum of A1 & A2 must be in the same plane
!        in order to keep the transverse momentum conservation.
        sctas2 = -ppt/ppa2
        cctas2 =  pz2/ppa2
        do i=1,3,1
            pa1(i) = pa(i)
            pa2(i) = pa(i)
        end do
        call rotate( cctas1, sctas1, cfis, sfis, ppa1, pa1, p00 )
!       Originally, momentum of A1 is relative to the direction of
!        the momentum of A
!       After rotation, momentum of A1 is relative to the CMS where A
!        belongs to.
        call rotate( cctas2, sctas2, cfis, sfis, ppa2, pa2, p00 )
!0130   Finished---------------------------------------

!0140   Newly induced parton (A2) should be added to particle list.
        nsca = nsca + 1
        kpip(nsca) = kf2
        pip(4,nsca) = amass(kf2)*amass(kf2)
        do i=1,3
            pip(4,nsca) = pip(4,nsca) + pa2(i)*pa2(i)
            pip(i,nsca) = pa2(i)
        end do
        pip(4,nsca) = SQRT( pip(4,nsca) )
!       pip(5,nsca) = qa2
        pip(5,nsca) = qqc
        pip(6,nsca) = xa2
!0140   Finished--------------------------------------------

!-------Updates five very important values of KFK, QA, XA, PA, EA (time-like)
!        & then perform the forward evolution along main chain further.
        kfk = kf1
        qa  = qqa
        do i=1,3
            pa(i) = pa1(i)
        end do
        if( i_split > 2 )then
            dm = amass(kf1)
        else
            dm = PSH(ij,5)
        end if
        ea = SQRT( pa(1)*pa(1) + pa(2)*pa(2) + pa(3)*pa(3) + dm*dm )
        xa = xa1

!       Records infomation of shifted mother and showered partons.
        ! Shifted mother or the first showerd parton.
        KSH( 2*ij+1, 1 ) = KSH(ij,1)
        KSH( 2*ij+1, 2 ) = kf1
        KSH( 2*ij+1, 3 ) = KSH(ij,3)
        KSH( 2*ij+1, 7 ) = KSH(ij,7)
        KSH( 2*ij+1, 8 ) = KSH(ij,8)
        ! 51: q -> qg; 52: g -> gg; 53: g -> qqbar
        KSH( 2*ij+1, 9 ) = i_split0 + 50
        ! 54: c -> cg; 55:        ; 56: g -> ccbar
        if( ABS( kf1 ) == 4 )then
            KSH( 2*ij+1, 9 ) = KSH( 2*ij+1, 9 ) + 3
        ! 57: b -> bg; 58:        ; 59: g -> bbbar
        else if( ABS( kf1 ) == 5 )then
            KSH( 2*ij+1, 9 ) = KSH( 2*ij+1, 9 ) + 6
        end if
        do i=1,3,1
            PSH( 2*ij+1, i ) = pa1(i)
        end do
        PSH( 2*ij+1, 4 ) = ea
        PSH( 2*ij+1, 5 ) = dm
        PSH( 2*ij+1, 6 ) = PSH(ij,6)
        PSH( 2*ij+1, 7 ) = PSH(ij,7)
        VSH( 2*ij+1, 1 ) = VSH(ij,1)
        VSH( 2*ij+1, 2 ) = VSH(ij,2)
        VSH( 2*ij+1, 3 ) = VSH(ij,3)
        VSH( 2*ij+1, 4 ) = VSH(ij,4)
        VSH( 2*ij+1, 5 ) = VSH(ij,4) + t_form1
        ! The second showered parton. Labels its status as "551".
        KSH( 2*ij+2, 1 ) = 551
        KSH( 2*ij+2, 2 ) = kf2
        KSH( 2*ij+2, 3 ) = KSH(ij,3)
        KSH( 2*ij+2, 9 ) = KSH( 2*ij+1, 9 )
        do i=1,3,1
            PSH( 2*ij+2, i ) = pa2(i)
        end do
        PSH( 2*ij+2, 4 ) = pip(4,nsca)
        PSH( 2*ij+2, 5 ) = amass(kf2)
        PSH( 2*ij+2, 6 ) = qqc
        PSH( 2*ij+2, 7 ) = 9D0
        VSH( 2*ij+2, 1 ) = VSH(ij,1)
        VSH( 2*ij+2, 2 ) = VSH(ij,2)
        VSH( 2*ij+2, 3 ) = VSH(ij,3)
        VSH( 2*ij+2, 4 ) = VSH(ij,4)
        VSH( 2*ij+2, 5 ) = VSH(ij,4) + t_form2

        i_succ = 1

!       if( n_sample_q > 10 ) return
        ! Once branching only.
        if( n_sample_q >= 1 ) return

!       i_stop=1 means the evolution "stop" while i_stop=0 "continue"
!       The time-like branching of A2 parton will consider later, i.e.
!        in "collis" after statement numbered 300.
        if( i_stop == 0 ) goto 10


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine suda_ti( i_fail, kf0, i_split, q2max, q2, z )
!!      Time-like Sudakov factor.
!!      Performs time-like forward branching, one step only.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (MSCA=20000)
        PARAMETER (pio=3.141592653589793D0)
        common/papr_p/core,xs,xu,xt,sm,as,dta,xa,sl0,tl0,qa, &
         ea,sqq,sqg,sgg,pa(3),pip(10,msca),mtime,kfk,nsca,kpip(msca)
        common/sa24/adj1(40),nnstop,non24,zstop
        real(kind=8) :: Lambda_QCD2
!----------------------------------------------------------------------
!       process    : A  ->   A1  +  A2
!       virtuality : qqb     qqa    qqc
!
!       equation: qqb = qqa / z + qqc / (1 - z) + pT2 / (1 - z) / z
!----------------------------------------------------------------------


!       yint1, 2: the integral functions of part of the splitting functions
        ! q -> q + g, cf. Sa' note p.44.
        yint1(x) = -LOG( 1D0 - x )
        ! g -> g + g
        yint2(x) =  LOG( 1D0 / (1D0 - x) - 1D0 )
        Lambda_QCD2 = adj1(25)*adj1(25)
        i_fail = 1
        i_split = 0
        q2 = 0D0
        z  = 0D0
        ! Cf. W' paper p.264 Eq.25.
        if( q2max <= 4*tl0 ) return
        ikf = ABS(kf0)
        tmax = LOG( q2max / Lambda_QCD2 )
        ! Cf. W' paper p.264
        tmin = LOG( 4D0*tl0 / Lambda_QCD2 )
!       tmax, tmin: bounds of t = LOG( Q^2/Lambda_QCD2 ) in Sudakov form factor
        ! cf. W' paper p.263
        zmin = tl0/q2max
        ! cf. W' paper p.263
        zmax = 1D0 - tl0/q2max
!       Approximated values of the bounds of z.
        if( zmax <= zmin ) return
        ! q -> qg
        if( ikf /= 21 )then
            ! cf. Sa' note p. 43
            ciqg = 4D0/3D0 * ( ( 2D0*yint1(zmax) - zmax - 0.5D0*zmax**2 ) &
                             - ( 2D0*yint1(zmin) - zmin - 0.5D0*zmin**2 ) )
!       ciqg: the integration of splitting function for q -> qg
            cisum = ciqg
            i_split = 1
        ! g ->
        else
            cigg = ( 6D0*yint2(zmax) - 12D0*zmax + 3D0*zmax**2 - 2D0*zmax**3 ) &
                 - ( 6D0*yint2(zmin) - 12D0*zmin + 3D0*zmin**2 - 2D0*zmin**3 )
!       cigg: the integration of the splitting function for g->gg
            ciqq = 1D0/2D0 * ( ( zmax - zmax**2 + 2D0/3D0*zmax**3 ) &
                             - ( zmin - zmin**2 + 2D0/3D0*zmin**3 ) )
!       ciqq: the integration of the splitting function for g -> qqbar
            cisum = cigg + ciqq
            aaa = PYR(1)
            ! g -> gg
            if( aaa <= ( cigg / cisum ) ) then
                i_split = 2
                cisum = cigg
            ! g -> qqbar
            else
                i_split = 3
                cisum = ciqq
            end if
        end if
        ce = 9D0/2D0/cisum

        n_sample_z = 0
100     continue

        aaa = PYR(1)
        ! Running alpha_s.
!       tt = tmax * aaa**ce
        ! Constant alpha_s.
        tt = tmax + 2D0*pio / as / cisum * LOG(aaa)
!       tt: the t value sampled from the time-like Sudakov factor
        if( tt <= tmin ) return
        q2 = Lambda_QCD2*EXP(tt)
!       Lambda_QCD2: Lambda_s**2
        rzmax = 0.5D0 + 0.5D0 * SQRT( 1D0 - 4D0*tl0/q2 )
!       cf. Eq.24 in B. R. Webber, Ann. Rev. Nucl. Part. Sci. 36(1986)253
!        that journal is simplified as W' elsewhere
        rzmin = 1D0 - rzmax
!       rzmax, rzmin: exact value of zmax & zmin corresponding to the t sampled.
!0170   Samples z. -------------------------------------------------------------
200     continue
        aaa = PYR(1)
        bbb = PYR(1)
        ! q -> qg
        if( i_split == 1 )then
            zsp = aaa*yint1(zmax) + (1D0 - aaa)*yint1(zmin)
            zsp = 1D0 - EXP(-zsp)
            ! The acceptance probability.
            ratio = (1D0 + zsp*zsp) / 2D0
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if( bbb > ratio ) goto 200
        ! g -> gg
        else if( i_split == 2 )then
            zsp = aaa*yint2(zmax) + (1D0 - aaa)*yint2(zmin)
            zsp = 1D0 - 1D0 / (1D0 + EXP(zsp) )
            ratio = ( 1D0 - zsp*(1D0 - zsp) )**2
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if( bbb > ratio ) goto 200
        ! g -> qqbar
        else
            zsp = aaa*zmax + (1D0 - aaa)*zmin
            ratio = 1D0 - 2D0*zsp + 2D0*zsp**2
            ! if(ratio <= bbb) goto 200
            ! Same as the previous statement.
            if(bbb > ratio) goto 200
        end if
!0170   Sampling z finished. ---------------------------------------------------

        if( zsp > rzmax .OR. zsp < rzmin ) then
!       If the z sampled falls out of the exact region, reject the sampled t.
!       Let tmax equal to the t sampled, and go back to resample.
            n_sample_z = n_sample_z + 1
!       Forcibly stops the time-like branching.
            if( n_sample_z > 1000 )then
                i_fail = 1
                return
            end if
            tmax = tt
            goto 100
        end if
        z = zsp
        i_fail = 0


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine cconse( ps, npl, np, pp, n_column, n_row )
!!      Adjusts four momentum conservation by iteration,no more than
!!       5000 iterations
!       pp : four momentum of particles
!       ps : above four momenta should be conserved to ps
!       npl : line # of the first particle
!       np : line # of last particle
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        real(kind=8) :: pp( n_column, n_row ), ff( n_column )
        dimension ps(4), pxyz(3), arp(3)


        ! Rescaling accuracy.
        ! dep  = 1D-5
        dep  = 1D-15
        ps4  = ps(4)
        pxyz = 0D0
        jj = 0
100     es = 0D0
        do i=npl,np,1
            if( pp(i,4) > 1D-10 ) es = es + pp(i,4)
        end do
        fr = es / ps4
!       if( ABS(1D0 - fr) <=  dep ) return
        n_rescale = 0
        do i=npl,np,1
            ppm   = pp(i,4)
            if( ppm <= 1D-10 ) cycle
            amas  = pp(i,5)
            amas2 = amas*amas
            ppf   = ppm / fr
            ff(i) = SQRT( ABS( ppf*ppf - amas2 ) / ( ppm*ppm - amas2 ) )
            do j=1,3,1
                ppp = ff(i) * pp(i,j)
                pp(i,j) = ppp
                pxyz(j) = pxyz(j) + ppp
            end do
            n_rescale = n_rescale + 1
        end do
        do i=1,3,1
            arp(i) = ABS( 1D0 - pxyz(i) / ps(i) )
        end do

        if( ABS(1D0 - fr) <= dep .AND. arp(1) <= dep .AND. arp(2) <= dep &
             .AND. arp(3) <= dep ) return

        do i=1,3,1
            pxyz(i) = pxyz(i) - ps(i)
            pxyz(i) = pxyz(i) / DBLE(n_rescale)
        end do
        do i=npl,np,1
            if( pp(i,4) <= 1D-10 ) cycle
            do j=1,3,1
                pp(i,j) = pp(i,j) - pxyz(j)
            end do
            pp5  = pp(i,5)
            pp52 = pp5 * pp5
            pp(i,4) = SQRT( pp52 + pp(i,1)**2 + pp(i,2)**2 + pp(i,3)**2 )
        end do
        jj = jj + 1
        if(jj < 5000) goto 100


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine perform_shower( icp )
!!      Performs the time-like shower (radiation), i.e. updates momenta and very
!!       important color tags of radiators, recoilers, and adds showerd paronts
!!       to the particle list /PYJETS/.
!       icp: "icp"-th entry in the time list.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        LOGICAL IS_QUARK, IS_DIQUARK
        PARAMETER (KSZJ=300000,KSZJ_PY8=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        COMMON/PYJET2/I_PTR2AA(KSZJ),K6K7K8(KSZJ,3),P6P7(KSZJ,2)
        COMMON/PYJET3/ N00_PY8, N00_DIQ, I_COL_MAX(2)
        common/shower_1/ KSH(14,10), PSH(14,10), VSH(14,10)
        common/shower_2/ NSH2, NSH2_max, IPTR_SH2(KSZJ), &
                         KSH2(KSZJ,8), PSH2(KSZJ,7), VSH2(KSZJ,5)
        common/sa6_p/ithroq_p,ithrob_p,ich_p,non6_p,throe_p(4)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/ctllist_p/nreac(100),nrel
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        common/work7/reac(100),crose(100)
!       Arrays of junction configuration of one NN pair for PYTHIA 8.
        INTEGER NO_JUNC_PY8, KIND_JUNC_PY8(KSZJ_PY8), &
                COL_JUNC_PY8(KSZJ_PY8,3), ENDC_JUNC_PY8(KSZJ_PY8,3), &
                STAT_JUNC_PY8(KSZJ_PY8,3)
        COMMON/PYJUNC_PY8/ NO_JUNC_PY8, KIND_JUNC_PY8, &
                           COL_JUNC_PY8, ENDC_JUNC_PY8, STAT_JUNC_PY8
!       lc(1,i), lc(2,i): line numbers in particle list of radiator and recoiler
!       lc(3,i): pointer to the current "NSH2"-th shower
!       lc(4,i): the shower case. =2, ic radiation; =3, jc radiation;
!                                 =12 and 13, ic & jc two radiations.
!       lc(5,i): the internal process number
!       lc(6,i): event type. =1, elastic scattering; =2, inelastic; =3, shower.
!       tc: the time of the shower will occur
!       ic, jc: line numbers of mothers in the particle list
!       ic_sh, jc_sh: line numbers of showered partons in the particle list
!       nreac(i): statistics of the # of successful i-th shower
        integer, intent(in) :: icp
        real(kind=8), dimension(4) :: ps
        real(kind=8), dimension(3,5) :: p_recoil

!       Location of strings.
        integer, parameter :: nn_str = 4000
        common/string_location/ n_string, i_type_string(nn_str), &
         i_string_PTR(KSZJ_PY8), n_in_string(nn_str), &
         i_location_str(1000,nn_str,2)


!       Gets contents of the current shower event.
        ic     = lc(1,icp)
        jc     = lc(2,icp)
        ic_sh  = 0
        jc_sh  = 0
        i_NSH2 = lc(3,icp)
        i_case = lc(4,icp)
        i_proc = lc(5,icp)
        i_SH2  = IPTR_SH2( i_NSH2 )

!       Adds the new radiated parton to "N+1". It is always stored in 2 of KSH2.
        i_mom = lc( i_case, icp )
        ! Note that 2 entries are a group for a shower event. Hence 2*.
        do i=1,5,1
            K( N+1, i ) = KSH2( 2 + 2*( i_SH2 - 1 ), i )
            P( N+1, i ) = PSH2( 2 + 2*( i_SH2 - 1 ), i )
            V( N+1, i ) = VSH2( 2 + 2*( i_SH2 - 1 ), i )
        end do
        K( N+1, 3 )      = i_mom - 1
        K6K7K8( N+1, 1 ) = 0
        P6P7( N+1, 1 )   = PSH2( 2 + 2*( i_SH2 - 1 ), 6 )
        P6P7( N+1, 2 )   = PSH2( 2 + 2*( i_SH2 - 1 ), 7 )
        V( N+1, 5 )      = 0D0

!       Sets color tags (K7 K8).
        KF_radiated = K( N+1, 2 )
        ! Gets the shower type.
        select case( i_proc )
        case( 51, 54, 57 )
            i_type_shower = 1
        case( 52 )
            i_type_shower = 2
        case( 53, 56, 59 )
            i_type_shower = 3
        case default
            i_type_shower = 0
            write(*,*) "Warning! perform_shower: unknown i_type_shower=", &
                       i_type_shower
            return
        end select
        ! Finds its alive ancestor (remaining), excepet the q broken from diq.
        N00   = N00_PY8
        i_anc = i_mom
        do while(.true.)
            KS_anc = K( i_anc, 1 )
            ! 543: a q/qbar broken from a (anti)diquark by PACIAE. Found.
            if( KS_anc == 543 )then
                i_anc = K( i_anc, 3 ) + 1
                exit
            else
                ! Found.
                if( KS_anc > 0 .AND. i_anc <= N00 ) exit
                ! Traces backward.
                i_anc = K( i_anc, 3 ) + 1
            end if
        end do
        KF_anc     = K( i_anc, 2 )
        i_col_anc  = K6K7K8( i_anc, 2 )
        i_acol_anc = K6K7K8( i_anc, 3 )
        ! New color/anti-color tag to be used.
        i_col_new = MAX( I_COL_MAX(1), I_COL_MAX(2) ) + 1

        ! Radiates a gluon: color and anti-color.
        if( KF_radiated == 21 )then
            ! q -> q + g
            if( i_type_shower == 1 )then
                ! Preares for the string connector (neighbor).
                i_connect  = 0
                KF_connect = 0
                i_col_connect  = -1
                i_acol_connect = -1
                ! If the ancestor is a (anti)diquark.
                if( IS_DIQUARK( KF_anc ) )then
                    ! If diquark, it has anti-color. Finds color.
                    if( KF_anc > 0 )then
                        i_col_diq = i_acol_anc
                        i_col = 2
                    ! If antidiquark, it has color. Finds anti-color.
                    else
                        i_col_diq = i_col_anc
                        i_col = 3
                    end if
                    ! Tries the backward searching at first.
                    do i = i_anc-1, 1, -1
                        KS_now = K(i,1)
                        if( KS_now <= 0 ) cycle
                        i_col_now = K6K7K8( i, i_col )
                        ! Finds its connector.
                        if( i_col_now == i_col_diq )then
                            i_connect  = i
                            KF_connect = K(i,2)
                            i_col_connect  = K6K7K8( i, 2 )
                            i_acol_connect = K6K7K8( i, 3 )
                            exit
                        end if
                    end do
                    ! If failed, tries the forward searching.
                    if( i_connect == 0 )then
                        ! do i = i_anc+1, N00, 1
                        do i = i_anc+1, N, 1
                            KS_now = K(i,1)
                            if( KS_now <= 0 ) cycle
                            i_col_now = K6K7K8( i, i_col )
                            ! Finds its connector.
                            if( i_col_now == i_col_diq )then
                                i_connect  = i
                                KF_connect = K(i,2)
                                i_col_connect  = K6K7K8( i, 2 )
                                i_acol_connect = K6K7K8( i, 3 )
                                exit
                            end if
                        end do
                    end if
                    ! Somthing goes wrong?
                    if( i_connect == 0 )then
                        write(*,*) "Warning! perform_shower-diq: " &
                                // "no color found, i_col_diq=", i_col_diq

                                !Lei_debug
                                write(22,*)
                                write(22,*) "Warning! perform_shower-diq: " &
                                         // "no color found, i_col_diq=", i_col_diq
                                write(22,*) "i_anc, KF_anc =", i_anc, KF_anc
                                write(22,*) "i_col_anc, i_acol_anc, i_col_new =", i_col_anc, i_acol_anc, i_col_new
                                write(22,*)
                                write(22,*) "i, type, col0, col1, col2, acol0, acol1, acol2:"
                                do ii_junc=1,NO_JUNC_PY8,1
                                     iijj_kind = KIND_JUNC_PY8( ii_junc )
                                     write(22,*) ii_junc-1, iijj_kind, &
                                         (  COL_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 ), &
                                         ( ENDC_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 )
                                end do
                                call PALIST(3)
                                write(22,*)
                                close(22)
                                stop
                                !Lei_debug
                        return
                    end if
                    ! Inserts this gluon to the diquark string.
                    if( i_col == 2 )then
                        K6K7K8( i_connect, 2 ) = i_col_new
                        K6K7K8( N+1, 2 ) = i_col_diq
                        K6K7K8( N+1, 3 ) = i_col_new
                    ! Anti-diquark string.
                    else
                        K6K7K8( i_connect, 3 ) = i_col_new
                        K6K7K8( N+1, 2 ) = i_col_new
                        K6K7K8( N+1, 3 ) = i_col_diq
                    end if

                ! If the ancestor is a (anti)quark.
                else if( IS_QUARK( KF_anc ) )then
                    if( KF_anc > 0 )then
                        i_col = 2
                    else
                        i_col = 3
                    end if
                    i_junc = 0
                    j_junc = 0
                    i_kind = 0
                    i_col_begin = -1
                    i_col_end   = -1
                    ! Number of junctions in the event.
                    N_JUNC = NO_JUNC_PY8
                    ! Ckecks if it is an end of a junction string.
                    loop_junc: do i = 1, N_JUNC, 1
                        do j=1,3,1
                            if( K6K7K8( i_anc, i_col ) &
                                == ENDC_JUNC_PY8(i,j)  )then
                                    i_junc = i
                                    j_junc = j
                                    i_kind = KIND_JUNC_PY8(i)
                                    i_col_begin = COL_JUNC_PY8(i,j)
                                    i_col_end   = ENDC_JUNC_PY8(i,j)
                                    exit loop_junc
                            end if
                        end do
                    end do loop_junc
                    ! If q, it has color. Finds anti-color.
                    if( KF_anc > 0 )then
                        i_col_q = i_col_anc
                        i_col = 3
                    ! If qbar, it has anti-color. Finds color.
                    else
                        i_col_q = i_acol_anc
                        i_col = 2
                    end if
                    i_col_junc = i_col_end

                    ! From a norma q-string.
                    if( i_junc == 0 )then
                        ! Tries the backward searching at first.
                        do i = i_anc-1, 1, -1
                            KS_now = K(i,1)
                            if( KS_now <= 0 ) cycle
                            i_col_now = K6K7K8( i, i_col )
                            ! Finds its connector.
                            if( i_col_now == i_col_q )then
                                i_connect  = i
                                KF_connect = K(i,2)
                                i_col_connect  = K6K7K8( i, 2 )
                                i_acol_connect = K6K7K8( i, 3 )
                                exit
                            end if
                        end do
                        ! If failed, tries the forward searching.
                        if( i_connect == 0 )then
                            ! do i = i_anc+1, N00, 1
                            do i = i_anc+1, N, 1
                                KS_now = K(i,1)
                                if( KS_now <= 0 ) cycle
                                i_col_now = K6K7K8( i, i_col )
                                ! Finds its connector.
                                if( i_col_now == i_col_q )then
                                    i_connect  = i
                                    KF_connect = K(i,2)
                                    i_col_connect  = K6K7K8( i, 2 )
                                    i_acol_connect = K6K7K8( i, 3 )
                                    exit
                                end if
                            end do
                        end if
                        ! Somthing goes wrong?
                        if( i_connect == 0 )then
                            write(*,*) "Warning! perform_shower-q-string:" &
                                    // " no color found, i_col_q=", i_col_q

                                    !Lei_debug
                                    write(22,*)
                                    write(22,*) "Warning! perform_shower-diq: " &
                                             // "no color found, i_col_q=", i_col_q
                                    write(22,*) "i_anc, KF_anc =", i_anc, KF_anc
                                    write(22,*) "i_col_anc, i_acol_anc, i_col_new =", i_col_anc, i_acol_anc, i_col_new
                                    write(22,*)
                                    write(22,*) "i, type, col0, col1, col2, acol0, acol1, acol2:"
                                    do ii_junc=1,NO_JUNC_PY8,1
                                         iijj_kind = KIND_JUNC_PY8( ii_junc )
                                         write(22,*) ii_junc-1, iijj_kind, &
                                             (  COL_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 ), &
                                             ( ENDC_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 )
                                    end do
                                    call PALIST(3)
                                    write(22,*)
                                    close(22)
                                    stop
                                    !Lei_debug
                            return
                        end if
                        ! Inserts this gluon to the q-string.
                        if( i_col == 2 )then
                            K6K7K8( i_connect, 2 ) = i_col_new
                            K6K7K8( N+1, 2 ) = i_col_q
                            K6K7K8( N+1, 3 ) = i_col_new
                        else
                            K6K7K8( i_connect, 3 ) = i_col_new
                            K6K7K8( N+1, 2 ) = i_col_new
                            K6K7K8( N+1, 3 ) = i_col_q
                        end if

                    ! From a junction string.
                    else
                        ! From a simple q(qbar) leg.
                        if( i_col_begin == i_col_end )then
                            ! Inserts this gluon directly and updates color tag.
                            if( i_col == 2 )then
                                K6K7K8( N+1, 2 ) = i_col_junc
                                K6K7K8( N+1, 3 ) = i_col_new
                            else
                                K6K7K8( N+1, 2 ) = i_col_new
                                K6K7K8( N+1, 3 ) = i_col_junc
                            end if
                            COL_JUNC_PY8( i_junc, j_junc ) = i_col_new
                        ! From a leg with gluons.
                        else
                            ! Tries the backward searching at first.
                            do i = i_anc-1, 1, -1
                                KS_now = K(i,1)
                                if( KS_now <= 0 ) cycle
                                i_col_now = K6K7K8( i, i_col )
                                ! Finds its connector.
                                if( i_col_now == i_col_junc )then
                                    i_connect  = i
                                    KF_connect = K(i,2)
                                    i_col_connect  = K6K7K8( i, 2 )
                                    i_acol_connect = K6K7K8( i, 3 )
                                    exit
                                end if
                            end do
                            ! If failed, tries the forward searching.
                            if( i_connect == 0 )then
                                ! do i = i_anc+1, N00, 1
                                do i = i_anc+1, N, 1
                                    KS_now = K(i,1)
                                    if( KS_now <= 0 ) cycle
                                    i_col_now = K6K7K8( i, i_col )
                                    ! Finds its connector.
                                    if( i_col_now == i_col_junc )then
                                        i_connect  = i
                                        KF_connect = K(i,2)
                                        i_col_connect  = K6K7K8( i, 2 )
                                        i_acol_connect = K6K7K8( i, 3 )
                                        exit
                                    end if
                                end do
                            end if
                            ! Somthing goes wrong?
                            if( i_connect == 0 )then
                                write(*,*) "Warning! perform_shower-j-string:" &
                                        // " no color found, i_col_junc=", &
                                           i_col_junc

                                           !Lei_debug
                                           write(22,*)
                                           write(22,*) "Warning! perform_shower-diq: " &
                                                    // "no color found, i_col_junc=", i_col_junc
                                           write(22,*) "i_anc, KF_anc =", i_anc, KF_anc
                                           write(22,*) "i_col_anc, i_acol_anc, i_col_new =", i_col_anc, i_acol_anc, i_col_new
                                           write(22,*)
                                           write(22,*) "i, type, col0, col1, col2, acol0, acol1, acol2:"
                                           do ii_junc=1,NO_JUNC_PY8,1
                                                iijj_kind = KIND_JUNC_PY8( ii_junc )
                                                write(22,*) ii_junc-1, iijj_kind, &
                                                    (  COL_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 ), &
                                                    ( ENDC_JUNC_PY8( ii_junc, jj_junc ), jj_junc=1,3,1 )
                                           end do
                                           write(22,*)
                                           call PALIST(3)
                                           write(22,*)
                                           close(22)
                                           stop
                                           !Lei_debug
                                return
                            end if
                            ! Inserts this gluon to the leg only.
                            if( i_col == 2 )then
                                K6K7K8( i_connect, 2 ) = i_col_new
                                K6K7K8( N+1, 2 ) = i_col_junc
                                K6K7K8( N+1, 3 ) = i_col_new
                            else
                                K6K7K8( i_connect, 3 ) = i_col_new
                                K6K7K8( N+1, 2 ) = i_col_new
                                K6K7K8( N+1, 3 ) = i_col_junc
                            end if
                        end if
                    end if

                end if

            ! g -> g + g
            else if( i_type_shower == 2 )then
                ! Inserts this gluon after the mother gluon directly.
                K6K7K8( i_anc, 3 ) = i_col_new
                K6K7K8( N+1, 2 ) = i_col_new
                K6K7K8( N+1, 3 ) = i_acol_anc
            else
                write(*,*) "Warning! perform_shower: " &
                        // "unkown i_type_shower-g =", i_type_shower
                return
            end if
        ! Quark, color.
        else if ( IS_QUARK( KF_radiated ) .AND. KF_radiated > 0 ) then
            !#TODO(Lei20250123): wating...
            return
        ! Anti-quark, anti-color.
        else if ( IS_QUARK( KF_radiated ) .AND. KF_radiated < 0 ) then
            !#TODO(Lei20250123): wating...
            return
        end if

!       Succeeded.

!       Sum of the initial 4-momenta.
        do i=1,4,1
            ps(i) = P(ic,i) + P(jc,i)
        end do

!       Updates momenta of the radiator. Adds a new radiated parton.
        select case( i_case )
        ! Case 1: "ic" radiates a parton "N+1".
        case(1)
            i_update = ic
            ic_sh = N + 1
            iic_sh = ic_sh
        ! Case 2: "jc" radiates a parton "N+1".
        case(2)
            i_update = jc
            jc_sh = N + 1
            iic_sh = jc_sh
        case default
            write(*,*) "Warning! perform_shower: unknown i_case1 =", i_case1
            return
        end select

        ! Updates momenta of the radiator "ic" or "jc". It is '1' in /shower_2/.
        do i=1,5,1
            P( i_update, i ) = PSH2( 1 + 2*( i_SH2 - 1 ), i )
        end do
        ! Adds the radiated parton correctly and enables it scat.
        N = N + 1
        ishp(N) = 1

!       Adjusts four momentum conservation by iterations, i.e. recoil.
        do j=1,5,1
            p_recoil( 1, j ) = P( ic, j )
            p_recoil( 2, j ) = P( jc, j )
            p_recoil( 3, j ) = P( iic_sh, j )
        end do
        n_column = 3
        n_row = 5
        call cconse( ps, 1, 3, p_recoil, n_column, n_row )
        ! Feeds the recoiled momenta back.
        do j=1,5,1
            P( ic, j )     = p_recoil( 1, j )
            P( jc, j )     = p_recoil( 2, j )
            P( iic_sh, j ) = p_recoil( 3, j )
        end do
        ! Collects potential lost 4-mom.
        do i=1,4,1
            throe_p(i) = throe_p(i) + ps(i) - P( ic, i )- P( jc, i ) &
                                            - P( iic_sh, i )
        end do

!       Updates flags, counters and pointers.
        reac(i_proc)  = reac(i_proc)  + 1D0
        crose(i_proc) = crose(i_proc) + 1D0
        nreac(i_proc) = nreac(i_proc) + 1
        I_COL_MAX = i_col_new
        IPTR_SH2( i_NSH2 ) = 0


        !Lei_debug
        ! Updates string locations
        i_str = i_string_PTR( i_anc )
        n_in_string( i_str ) = n_in_string( i_str ) + 1
        n_in_str = n_in_string( i_str )
        i_location_str( n_in_str, i_str, 2 ) = N
        !Lei_debug


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine his_p( ic, jc, ic_sh, jc_sh, pi00, pj00, time_in, i_type )
!!      Propagates particles to the collision/shower time in the Lab frame.
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        integer, intent(in) :: ic, jc, ic_sh, jc_sh, i_type
        real(kind=8), intent(in) :: pi00(4), pj00(4), time_in
        ! pi00, pj00: four-momenta before the collisions/shower
        ! i_type: event type. =1, elastic scattering; =2, inelastic; =3, shower.


        !Lei_debug
        ! Energy loss and Mean free path.
        common/mean_path/ path_sum(KSZJ), coll_sum(KSZJ)
        real(kind=8) :: vi0(3), vj0(3), dxloss(2)
        do i=1,3,1
            vi0(i) = V(ic,i)
            vj0(i) = V(jc,i)
        end do
        !Lei_debug


!       Propagates scattered particles or the radiator and the radiated parton.
        ! Scattering event.
        if( i_type == 1 .OR. i_type == 2 )then
            ! ic
            tau_p = 0D0
            ! Tries to give a formation time for the inelastic scattering (?).
            if( i_type == 2 ) tau_p = 0D0
            time = time_in + tau_p
            do j=1,3,1
                V(ic,j) = V(ic,j) + pi00(j) / pi00(4) * ( time - V(ic,4) )
            end do
            V(ic,4)  = time
            taup(ic) = V(ic,4)
            ! jc
            tau_p = 0D0
            ! Tries to give a formation time for the inelastic scattering (?).
            if( i_type == 2 ) tau_p = 0D0
            time = time_in + tau_p
            do j=1,3,1
                V(jc,j) = V(jc,j) + pj00(j) / pj00(4) * ( time - V(jc,4) )
            end do
            V(jc,4)  = time
            taup(jc) = V(jc,4)
        ! Shower event.
        else if( i_type == 3 )then
            ! Paton shower from "ic".
            if( ic_sh > 0 )then
                ! Propagates the radiator.
                time = time_in
                do j=1,3,1
                    V(ic,j) = V(ic,j) + pi00(j) / pi00(4) * ( time - V(ic,4) )
                end do
                V(ic,4)  = time
                taup(ic) = V(ic,4)
                ! Propagates the radiated parton.
                do j=1,4,1
                    V( ic_sh, j ) = V( ic, j )
                end do
                taup( ic_sh ) = V( ic_sh, 4 )
            end if
            ! Paton shower from "jc".
            if( jc_sh > 0 )then
                ! Propagates the radiator.
                time = time_in
                do j=1,3,1
                    V(jc,j) = V(jc,j) + pj00(j) / pj00(4) * ( time - V(jc,4) )
                end do
                V(jc,4)  = time
                taup(jc) = V(jc,4)
                ! Propagates the radiated parton.
                do j=1,4,1
                    V( jc_sh, j ) = V( jc, j )
                end do
                taup( jc_sh ) = V( jc_sh, 4 )
            end if
        end if


        !Lei_debug
        dxloss(1) = SQRT( pow2( vi0(1) - V(ic,1) ) + pow2( vi0(2) - V(ic,2) ) + pow2( vi0(3) - V(ic,3) ) )
        dxloss(2) = SQRT( pow2( vj0(1) - V(jc,1) ) + pow2( vj0(2) - V(jc,2) ) + pow2( vj0(3) - V(jc,3) ) )
        path_sum(ic) = path_sum(ic) + dxloss(1)
        path_sum(jc) = path_sum(jc) + dxloss(2)
        coll_sum(ic) = coll_sum(ic) + 1D0
        coll_sum(jc) = coll_sum(jc) + 1D0
        !Lei_debug


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine update_ctl( icp )
!!      Throws away old collission pairs and updates the time list.
!       Imposes the LPM effect of the shower event via the time constraint.
!       icp: scattered "icp"-th collision pair
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/shower_1/ KSH(14,10), PSH(14,10), VSH(14,10)
        common/shower_2/ NSH2, NSH2_max, IPTR_SH2(KSZJ), &
                         KSH2(KSZJ,8), PSH2(KSZJ,7), VSH2(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        integer, intent(in) :: icp
        integer :: i_NSH2(2)
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


!       Throws away the current scattering pair.
        icol = icol - 1
        do k1 = icp, icol, 1
            lc( 1, k1 ) = lc( 1, k1+1 )
            lc( 2, k1 ) = lc( 2, k1+1 )
            lc( 3, k1 ) = lc( 3, k1+1 )
            lc( 4, k1 ) = lc( 4, k1+1 )
            lc( 5, k1 ) = lc( 5, k1+1 )
            lc( 6, k1 ) = lc( 6, k1+1 )
            tc( 1, k1 ) = tc( 1, k1+1 )
            tc( 2, k1 ) = tc( 2, k1+1 )
            sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
            sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
            sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
            sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
        end do

!       Records the information of showers in /shower_1/ into /shower_2/ if any.
        i_NSH2 = -999
        i_shower_survived = 0
        ! Up to max two showers only!
        iNSH = NSH2
        do i_shower = 1, 2, 1
            ! Position of 1st shower in /shower_1/.
            iSH = 4
            ! 2nd
            if( i_shower == 2 ) iSH = 6
            ! 551 : shower status of PACIAE parcas.
            if( KSH( iSH, 1 ) /= 551 ) cycle
            ! Records the radiator "iSH-1" and the radiated parton "iSH" to
            !  the "iNSH + 1"-th group of entries in /shower_2/.
            iNSH = iNSH + 1
            i_NSH2( i_shower ) = iNSH
            ii = 0
            ! NB: 2 entries as a group. So 2*.
            do i = iSH-1, iSH, 1
                ii = ii + 1
                do j=1,5,1
                    KSH2( ii + 2*( iNSH - 1 ), j ) = KSH(i,j)
                    PSH2( ii + 2*( iNSH - 1 ), j ) = PSH(i,j)
                    VSH2( ii + 2*( iNSH - 1 ), j ) = VSH(i,j)
                end do
                do j=6,7,1
                    KSH2( ii + 2*( iNSH - 1), j ) = KSH(i,j)
                    PSH2( ii + 2*( iNSH - 1), j ) = PSH(i,j)
                end do
                ! KSH2( ii + 2*iNSH, 8 ) = KSH(i,8)
                KSH2( ii + 2*( iNSH - 1 ), 8 ) = KSH(i,8)
            end do
        end do

!       Loops over ic, jc (new) and old partons (i.e. constructs colli. pair
!        by partons, one of which is ic/jc and another one is in parton list).

        ! Time resolution.
        dddt  = adj1(19)
        t_max = adj1(28)
        ! The collision pair counter (+1 avoids 0 icol).
        icol = icol + 1
        do ii=1,2,1
            if( ii == 1 ) i = ic
            if( ii == 2 ) i = jc
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = 1, N_old, 1
                if( ishp(j) == 0 ) cycle
                ! Avoids two particles colliding immediately after a collision.
                if( j == ic .OR. j == jc ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    i_case = lc(4,j1)
                    i_type = lc(6,j1)
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    ! Shower event of "j1".
                    if( i_type == 3 )then
                        if( i_case == 1 )then
                            jjc = iic
                        elseif( i_case == 2 )then
                            iic = jjc
                        ! Should not happen.
                        else
                            write(*,*) "Warning! update_ctl: unknow " &
                                    // "i_case-1:", i_case
                            iic = 0
                            jjc = 0
                        end if
                    end if
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Throws away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = 2
                        ! For the shower event, clears the pointers too.
                        if( i_type == 3 )then
                            ii_NSH2 = lc(3,j1)
                            IPTR_SH2( ii_NSH2 ) = 0
                        end if
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time calculation.
                icol = icol + 1
                ! Tries to add a shower event.
                if( i_NSH2( ii ) < 0 ) cycle
                ! Checks if "icol" is the one related to 'i'.
                ! 'i' related pair has been thrown away. The shower might occur.
                iic = lc( 1, icol-1 )
                ttc = tc( 1, icol-1 )
                if( iic /= i )then
                    t_coll = t_max
                else
                    t_coll = ttc
                end if
                if( i == ic )then
                    iSH = 4
                else
                    iSH = 6
                end if
                ! Gets the shower type.
                i_proc = KSH( iSH, 9 )
                select case( i_proc )
                case( 51, 54, 57 )
                    i_type_shower = 1
                case( 52 )
                    i_type_shower = 2
                case( 53, 56, 59 )
                    i_type_shower = 3
                case default
                    i_type_shower = 0
                    write(*,*) "Warning! update_ctl: unknown i_type_shower-1=",&
                               i_type_shower
                    cycle
                end select
                ! With potential showers of "i". Only q -> q + g and g -> g + g.
                if( i_type_shower == 1 .OR. i_type_shower == 2 )then
                    t_shower = VSH( iSH, 5 )
                    ! Scattering will suppress (kill) this shower event. LPM.
                    if( t_shower >= t_coll )then
                        ! Do nothing.
                    ! Survives. Replaces the "icol-1"-th scattering event of i"
                    !  with this shower event in the time list or adds new one.
                    else
                        if( iic /= i )then
                            k1 = icol
                        else
                            k1 = icol - 1
                        end if
                        lc( 1, k1 ) = ic
                        lc( 2, k1 ) = jc
                        ! Now "lc3" points to "NSH2".
                        NSH2 = i_NSH2(ii)
                        lc( 3, k1 ) = NSH2
                        IPTR_SH2( NSH2 ) = NSH2
                        ! Now "lc4" labels the shower case.
                        if( i == ic )then
                            lc( 4, k1 ) = 1
                        else
                            lc( 4, k1 ) = 2
                        end if
                        lc( 5, k1 ) = i_proc
                        lc( 6, k1 ) = 3
                        tc( 1, k1 ) = t_shower
                        tc( 2, k1 ) = t_shower
                        i_shower_survived = i_shower_survived + lc( 4, k1 )
                        ! For the next time calculation.
                        if( iic /= i )then
                            icol = icol + 1
                        end if
                    end if
                end if
            ! No scattering event of 'i' will occur.
            ! Adds a shower event into the time list if any.
            else
                if( i_NSH2( ii ) < 0 ) cycle
                if( i == ic )then
                    iSH = 4
                else
                    iSH = 6
                end if
                ! Gets the shower type.
                i_proc = KSH( iSH, 9 )
                select case( i_proc )
                case( 51, 54, 57 )
                    i_type_shower = 1
                case( 52 )
                    i_type_shower = 2
                case( 53, 56, 59 )
                    i_type_shower = 3
                case default
                    i_type_shower = 0
                    write(*,*) "Warning! update_ctl: unknown i_type_shower-2=",&
                               i_type_shower
                    cycle
                end select
                ! With potential showers of "i". Only q -> q + g and g -> g + g.
                if( i_type_shower == 1 .OR. i_type_shower == 2 )then
                    t_shower = VSH( iSH, 5 )
                    if( t_shower > t_max )then
                        ! Do nothing.
                    else
                        k1 = icol
                        lc( 1, k1 ) = ic
                        lc( 2, k1 ) = jc
                        ! Now "lc3" points to "NSH2".
                        NSH2 = i_NSH2(ii)
                        lc( 3, k1 ) = NSH2
                        IPTR_SH2( NSH2 ) = NSH2
                        ! Now "lc4" labels the shower case.
                        if( i == ic )then
                            lc( 4, k1 ) = 1
                        else
                            lc( 4, k1 ) = 2
                        end if
                        lc( 5, k1 ) = i_proc
                        lc( 6, k1 ) = 3
                        tc( 1, k1 ) = t_shower
                        tc( 2, k1 ) = t_shower
                        i_shower_survived = i_shower_survived + lc( 4, k1 )
                        ! For the next time calculation.
                        icol = icol + 1
                    end if
                end if
            end if
        end do
        ! Deducts 1 because the counter begins with +1 not 0.
        icol  = icol - 1
        N_old = N


        return
        end



!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        subroutine update_time_list_shower( icp )
!!      Throws away old collission and shower events. Updates the time list.
!       Imposes the LPM effect of the shower event via the time constraint.
!       icp: showered "icp"-th event
        IMPLICIT DOUBLE PRECISION(A-H, O-Z)
        IMPLICIT INTEGER(I-N)
        PARAMETER (KSZJ=300000,MCLIS=280000)
        COMMON/PYJETS/N,NPAD,K(KSZJ,5),P(KSZJ,5),V(KSZJ,5)
        common/shower_1/ KSH(14,10), PSH(14,10), VSH(14,10)
        common/shower_2/ NSH2, NSH2_max, IPTR_SH2(KSZJ), &
                         KSH2(KSZJ,8), PSH2(KSZJ,7), VSH2(KSZJ,5)
        common/sa8_p/taup(kszj),coor(3),ishp(kszj)
        common/sa24/adj1(40),nnstop,non24,zstop
        common/collist/lc(6,mclis),tc(2,mclis),sig_tsmp_Cthetas_p(4,mclis),icol
        common/scatt/ ic, jc, ic_sh, jc_sh, N_old
        integer, intent(in) :: icp
        integer :: l_coll(6)
        real(kind=8) :: t_coll, t_coll_min
        real(kind=8) :: sig_t_theta_p(4)
        logical :: succeed


!       Throws away the currently shower event.
        icol = icol - 1
        do k1 = icp, icol, 1
            lc( 1, k1 ) = lc( 1, k1+1 )
            lc( 2, k1 ) = lc( 2, k1+1 )
            lc( 3, k1 ) = lc( 3, k1+1 )
            lc( 4, k1 ) = lc( 4, k1+1 )
            lc( 5, k1 ) = lc( 5, k1+1 )
            lc( 6, k1 ) = lc( 6, k1+1 )
            tc( 1, k1 ) = tc( 1, k1+1 )
            tc( 2, k1 ) = tc( 2, k1+1 )
            sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
            sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
            sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
            sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
        end do

!       Loops over ic, jc and ic_sh/jc_sh (new) and old partons (i.e. constructs
!        colli. pair by partons, one of which is ic/jc/ic_sh/jc_sh and another
!        one is in the old parton list).

        ! Time resolution.
        dddt  = adj1(19)
        t_max = adj1(28)
        ! The collision pair counter (+1 avoids 0 icol).
        icol = icol + 1
        do ii=1,3,1
            if( ii == 1 ) i = ic
            if( ii == 2 ) i = jc
            ! The showerd parton should be in 'N'.
            if( ii == 3 ) i = N
            if( ishp(i) == 0 ) cycle
            ! Successful flag.
            succeed = .false.
            ! A number large enough to determine the minumum time.
            t_coll_min = 1D30
            ixc = 0
            jxc = 0
            loop_j: do j = 1, N_old, 1
                if( ishp(j) == 0 ) cycle
                ! Avoids particles colliding immediately after shower.
                if( j == ic .OR. j == jc ) cycle
                call coij_p( i, j, i_fail, l_coll, t_coll, sig_t_theta_p )
                if( i_fail == 1 ) cycle
                ! Imposes the time resolution constraint.
                if( t_coll < 1D-10 ) cycle
                do j1 = 1, icol-1, 1
                    if( ABS( tc(1,j1) - t_coll ) < dddt ) cycle loop_j
                end do
                ! Chooses the smallest time for 'i'-cycle.
                if( t_coll < t_coll_min )then
                    lc(1,icol) = l_coll(1)
                    lc(2,icol) = l_coll(2)
                    lc(3,icol) = l_coll(3)
                    lc(4,icol) = l_coll(4)
                    lc(5,icol) = l_coll(5)
                    lc(6,icol) = l_coll(6)
                    tc(1,icol) = t_coll
                    tc(2,icol) = t_coll
                    sig_tsmp_Cthetas_p(1,icol) = sig_t_theta_p(1)
                    sig_tsmp_Cthetas_p(2,icol) = sig_t_theta_p(2)
                    sig_tsmp_Cthetas_p(3,icol) = sig_t_theta_p(3)
                    sig_tsmp_Cthetas_p(4,icol) = sig_t_theta_p(4)
                    t_coll_min = t_coll
                    ixc = l_coll(1)
                    jxc = l_coll(2)
                end if
                succeed = .true.
            end do loop_j
            ! Keeps the one with the smallest time from pairs including i or j.
            if( succeed )then
                n_jump_out = 0
                j1 = 1
                do while(.true.)
                    if( j1 > icol-1 ) exit
                    i_case = lc(4,j1)
                    i_type = lc(6,j1)
                    iic = lc(1,j1)
                    jjc = lc(2,j1)
                    ! Shower event of "j1".
                    if( i_type == 3 )then
                        if( i_case == 1 )then
                            jjc = iic
                        elseif( i_case == 2 )then
                            iic = jjc
                        ! Should not happen.
                        else
                            write(*,*) "Warning! update_time_list_shower: " &
                                    // "unknow  i_case:", i_case
                            iic = 0
                            jjc = 0
                        end if
                    end if
                    if(       ixc /= iic .AND. ixc /= jjc &
                        .AND. jxc /= iic .AND. jxc /= jjc )then
                        j1  = j1 + 1
                        cycle
                    end if
                    ttc = tc(1,j1)
                    ! Throws away the pair with larger time.
                    if( ttc > t_coll_min )then
                        k_begin = j1
                        n_jump_out = 2
                        ! For the shower event, clears the pointers too.
                        if( i_type == 3 )then
                            ii_NSH2 = lc(3,j1)
                            IPTR_SH2( ii_NSH2 ) = 0
                        end if
                    else
                        k_begin = icol
                        n_jump_out = 2
                    end if
                    icol = icol - 1
                    do k1 = k_begin, icol, 1
                        lc( 1, k1 ) = lc( 1, k1+1 )
                        lc( 2, k1 ) = lc( 2, k1+1 )
                        lc( 3, k1 ) = lc( 3, k1+1 )
                        lc( 4, k1 ) = lc( 4, k1+1 )
                        lc( 5, k1 ) = lc( 5, k1+1 )
                        lc( 6, k1 ) = lc( 6, k1+1 )
                        tc( 1, k1 ) = tc( 1, k1+1 )
                        tc( 2, k1 ) = tc( 2, k1+1 )
                        sig_tsmp_Cthetas_p(1,k1) = sig_tsmp_Cthetas_p( 1, k1+1 )
                        sig_tsmp_Cthetas_p(2,k1) = sig_tsmp_Cthetas_p( 2, k1+1 )
                        sig_tsmp_Cthetas_p(3,k1) = sig_tsmp_Cthetas_p( 3, k1+1 )
                        sig_tsmp_Cthetas_p(4,k1) = sig_tsmp_Cthetas_p( 4, k1+1 )
                    end do
                    if( n_jump_out == 2 ) exit
                end do
                ! For the next time calculation.
                icol = icol + 1
            end if
        end do
        ! Deducts 1 because the counter begins with +1 not 0.
        icol  = icol - 1
        N_old = N


        return
        end



!ccccccccccccccccccccccccccccccccccccc end ccccccccccccccccccccccccccccccccccccc